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The  study  reported  in  this  thesis  deals  with  the 
development  of  an  automatic  shape  optimization  procedure  of 
two-dimensional  structures  for  any  given  loads,  support 
conditions  and  geometric  constraints.  The  procedure  utilizes 
a specially  developed  finite  element  program  with  the 
appropriate  automated  rule-based  mesh  generation  for  the 
structural  modeling  and  stress  evaluation  at  each  step  of  the 
optimization  process. 

The  process  starts  with  a primitive  shape  and  continues 
to  change  it  by  eliminating  groups  of  elements  at  each  step 
according  to  the  preprogrammed  rules.  The  final  shape  is 
determined  when  no  further  elimination  is  possible  without 
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violating  the  allowable  stress,  or  when  the  average  stress 
reaches  the  target  value  relative  to  maximum  stress  in  the 
entire  structure. 

The  developed  procedure  is  illustrated  by  application  to 
several  examples  of  2-D  structures  subjected  to  static  loads, 
inertial  (body)  forces  and  impact.  The  results  demonstrate 
the  excellent  potential  of  the  approach  for  general  use  in 
automated  shape  optimization  under  any  given  loads,  boxindary 
conditions  and  imposed  constraints. 


CHAPTER  1 
INTRODUCTION 

In  mechanical  design,  there  are  a large  number  of 
applications  in  which  the  engineer's  primary  concern  is  the 
optimum  utilization  of  materials.  The  automotive  and  aircraft 
industries  applied  this  concept  to  many  part  designs  in  order 
to  reduce  weight  and  conseguently  conserve  the  energy  used  by 
the  vehicles.  There  are  many  ways  to  approach  the  target  of 
fully  stressed  design  in  the  field  of  shape  optimization  of 
mechanical  components  or  structural  elements.  The  purpose  of 
all  these  considerations  is  to  produce  a design  in  which  each 
structural  element  sustains  a stress  approaching  the  allowable 
stress  of  material  under  the  given  loads  and  support 
conditions. 

Problems  of  shape  optimization  arise  in  situations  where 
the  shape  parameters  are  computed  as  unknowns  during  the 
optimization  process.  Selection  of  design  variables  is  case 
dependent  and  the  result  of  the  optimization  procedures  may 
not  be  unique.  As  the  structure  and  the  boundary  conditions 
become  complex,  it  is  not  easy  to  develop  a general  algorithm 
that  provides  an  optimum  configuration.  Furthermore,  the 
computational  time  needed  in  the  optimization  process 
increases  rapidly  as  the  number  of  the  design  variables 
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increases.  There  are  many  pxiblished  methods  to  solve  the 
shape  optimization  problem,  and  we  will  briefly  review  some  of 
them  in  the  following. 

Different  Approaches  of  Form  Synthesis 

In  form  synthesis  the  optimum  structure  can  be  attained 
in  two  ways.  The  first  is  based  on  changing  the  material 
properties  such  as  elastic  modulus,  density  and  strength 
without  changing  the  shape  of  the  body.  The  second  is  based 
on  changing  the  shape  without  changing  the  material 
properties.  This  study  belongs  to  the  second  classification. 
The  shape  of  the  structure  is  continuously  changed  in  the 
design  process  until  shape  optimization  with  minimum  weight  is 
achieved.  Shape  optimization  of  structural  elements  or 
mechanical  components,  in  addition  to  weight  reduction,  is 
frequently  incorporated  with  secondary  objectives  such  as 
minimization  of  stress  concentration,  reduction  of  maximum 
stress  and  improvement  of  fatigue  behavior. 

Mathematical  Programming 

In  the  early  1960s  the  numerical  approximation  techniques 
were  combined  with  mathematical  programming  methods  to  develop 
structural  optimization  algorithms.  The  advances  in 
mathematical  programming  techniques  and  high  speed  computers 
gave  an  impetus  to  the  wide  application  of  optimization  in 
engineering  design. 

When  dealing  with  mathematical  programming  procedures  in 
complex  problems,  powerful  tools  are  usually  needed  to 
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determine  the  magnitudes  of  displacement,  stress,  force  or 
energy  that  are  involved  in  the  calculation  of  the  objective 
function.  The  discretization  modeling  procedures  that  are 
frequently  used  in  this  kind  of  approach  are  the  finite 
element  method  and  the  boundary  element  method. 

Applications  of  programming  and  variational  techniques  in 
the  mechanical  systems  have  been  mainly  in  the  fields  of 
structural  design,  synthesis  of  control  systems  and 
optimization  of  rocket  and  space  vehicle  performance.  Seireg 
[1]  presented  a review  of  some  illustrative  examples  of  the 
uses  of  optimization  techniques  in  the  design  of  mechanical 
elements  and  systems.  The  cases  cited  included  gears, 
bearings,  dynamic  systems,  rotating  discs,  pressure  vessels, 
shaft  under  bending  and  torsion,  beams  subjected  to 
longitudinal  impact,  and  problems  of  elastic  contact  and  load 
distribution. 

Belegundu  and  Arora  [2]  presented  a comprehensive  review 
of  various  mathematical  programming  techniques  used  for 
structural  optimization.  The  applicability  of  modern 
optimization  techniques  to  structural  design  problems  was 
discussed  from  the  design  engineers'  viewpoint. 

Although  many  optimization  problems  in  engineering  design 
are  not  strictly  linear,  the  engineer  practically  always  uses 
approximate  optimization  based  on  simple  linearized  model 
rather  than  the  exact  solution,  since  linear  problems  are 
easer  to  solve  than  nonlinear  ones.  R.  K.  Livesley  [3] 
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introduced  the  linear  programming  method  in  structural 
analysis  and  design.  The  most  prevalent  applications  of 
linear  programming  in  structural  engineering  have  occurred  in 
limit  analysis  and  design,  where  the  designer  looks  for  upper- 
and  lower-bounds  solutions. 

In  many  optimization  processes  the  determination  of 
sensitivities  of  the  objective  function  and  the  constraints  to 
change  in  the  design  variables  are  important  factors  in  the 
procedures.  In  general,  only  certain  explicit  constraints  are 
simple  linear  functions  of  the  design  variables,  and  both  the 
objective  function  and  implicit  constraints  are  dependent  in 
a nonlinear  manner  on  the  design  variables.  Zienkiewicz  and 
Campbell  [4]  devised  an  efficient  method  of  determining  the 
sensitivities  of  stress/deformation  guantities  under  some 
circumstances  using  a sequence  of  linear-programming 
operations . 

Fleury  [5]  presented  the  convex  linearization  method  for 
shape  optimal  design  of  elastic  structure  discretized  in 
finite  element.  This  approach  is  based  on  an  internal 
parametric  representation  of  modern  techniques  employed  in 
computer  aided  geometric  design.  The  designed  geometric 
boundary  is  controlled  by  master  nodes.  This  method  makes  the 
explicit  constraints  convex  by  partially  linearizing  them  with 
respect  to  the  reciprocal  variables.  The  objective  function 
was  linearized  with  respect  to  the  direct  variables  and  the 
constraints  were  linearized  with  respect  to  the  reciprocal 
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When  the  objective  function  and  constraints  are 
nonlinear,  as  in  most  cases  of  structural  design,  it  is  not 
possible  to  use  the  linear  programming  methods.  In  practical 
structural  design  problems,  the  solutions  are  usually  located 
on  the  border  of  the  feasible  region  by  the  imposed  stress  and 
deflection  constraints.  This  characteristic  of  the  optimum 
design  is  utilized  in  many  optimization  algorithms  that 
involve  systematic  searches  for  the  best  solution  at  or  in  the 
vicinity  of  the  boundary  between  the  feasible  and  the 
infeasible  regions.  Queau  and  Trompette  [6]  used  an  extended 
interior  penalty  function  to  find  the  optimal  shape  of  plane 
or  axisymmetric  structures  in  order  to  minimize  the  stress 
concentration  factor  along  the  boundary  described  by  straight 
lines  and  circles.  A general  process,  which  took  into  account 
some  developments  in  the  field  of  numerically  controlled 
machines , was  presented  in  which  the  entire  boundary  can  be 
described  without  any  constraints  on  the  change  scale. 

Dems  and  Mroz  [7]  presented  a multiparameter  structural 
shape  optimization  using  the  finite  element  method  and  an 
optimality  criterion.  The  problem  of  optimal  design  of  the 
shape  of  a free  or  internal  boundary  of  a body  is  formulated 
by  assuming  that  the  boundary  shape  is  described  by  a set  of 
prescribed  shape  functions  and  a set  of  parameters.  The 
optimization  procedure  is  reduced  to  the  determination  of 
these  parameters.  This  allows  the  change  of  boundary 
conditions  with  the  shape  modification  that  is  needed  in  the 
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problem  of  determination  of  shape  of  a free  or  loaded  boundary 
of  a plate  or  configuration  of  truss  or  frame  structure.  The 
optimality  criteria  approach  essentially  involves  the 
derivation  of  necessary  and  sufficient  conditions  for  an 
optimum  solution  with  specific  design  requirements, 
associating  them  with  a certain  real  or  pseudo  energy 
functional  of  the  system.  The  finite  element  method  naturally 
determines  many  of  these  energy  functions,  and  there  is  no 
additional  transformation  needed  for  determining  the 
directions  in  the  searching  process. 

Geometric  Approaches 

With  computer-aided  geometric  design  systems  becoming 
more  and  more  popular  in  design  applications,  it  is  a natural 
thing  to  couple  the  CAD  system  with  structural  optimization. 
Esping  [8]  used  a CAD  approach  to  the  minimum  weight  design 
problem.  A structure  is  defined  by  a set  of  governing  points. 
These  points  can  be  used  to  define  lines,  surfaces  and  solids. 
A governing  point  can  be  connected  to  design  variables.  The 
values  of  the  design  variables  are  determined  in  order  to 
minimize  the  weight  of  the  structure  subject  to  a number  of 
constraints  on  displacements  or  stresses . A truss  structure 
can  be  described  as  a set  of  lines,  a membrane  or  shell 
structure  as  a set  of  surfaces,  and  solid  structures  as 
solids.  In  all  cases  points  are  the  basic  information.  A new 
shape  is  created  by  moving  the  points.  It  is  then  natural  to 
connect  design  variables  to  the  point  coordinates. 
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Braibant  and  Fleury  [9]  also  used  the  CAD  philosophy  of 
geometric  description  to  formulate  shape  optimal  design  of  an 
elastic  structure.  The  structure  to  be  optimized  is 
decomposed  into  a set  of  simple  subregions,  which  are  called 
design  elements,  and  the  shape  of  these  design  elements  is 
described  by  some  properly  chosen  master  nodes.  Parametric 
curves  and  surfaces  of  CAD  realize  the  interpolation  process 
in  the  design  elements.  The  master  nodes  positions  are  the 
design  variables.  The  advantage  of  this  approach  is  that  it 
is  no  longer  necessary  either  to  piece  design  elements 
together  so  as  to  agree  with  the  shape  complexity  or  to 
restrict  the  shape  variations.  Wang,  Sun  and  Gallagher  [10] 
also  used  this  kind  of  geometric  description  in  the 
sensitivity  analysis  for  shape  optimization  of  continuum 
structure . 

Bennett  and  Botkin  [11]  used  automatic  mesh  generation 
that  was  based  only  on  boundary  points  to  couple  with  adaptive 
refinement  in  the  shape  optimization  of  plate  structures.  The 
structure  is  modeled  by  using  boundary  design  elements  and  key 
nodes.  The  geometry  of  the  elements  is  defined  in  two 
dimensions  by  shape  functions  along  with  design  variables. 
The  boundary  elements  are  also  used  to  define  the  stress 
constraints.  Each  boundary  design  element  is  associated  with 
at  least  one  stress  constraint  which  is  calculated  from  the 
maximum  stress  of  all  finite  elements  touching  that  boundary 
design  element.  The  loads  and  structural  boundary  conditions 
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are  related  to  a set  of  reference  nodes,  which  are  called  key 
nodes.  The  shape  of  the  structure  is  changed  by  adaptation  of 
the  shape  of  the  boundary  design  elements. 

A study  by  Imam  [12]  deals  directly  with  three- 
dimensional  shape  optimization  of  structural  components  by  the 
concept  of  shape  algorithm.  In  the  context  of  shape 
optimization,  the  objective  is  to  change  the  shape  of  a 
surface  by  changing  parameters  called  the  shape  variables  or 
the  design  variables.  If  some  parametric  variables  can 
identify  the  points  on  a surface,  then  the  numerical 
representation  of  a surface  requires  a method  to  compute  the 
spatial  coordinates  of  any  point  on  the  surfaces  from  the 
shape  variables  and  parametric  variables.  Mass  is  minimized 
by  using  the  feasible  direction  method.  Numerical  shape 
representation  and  the  selection  of  design  variables  are 
important  aspects  treated  by  means  of  isoparametric 
representation  of  the  surface  and  by  numerical  superposition 
of  shapes.  In  general,  the  problem  is  addressed  from  a 
practical  standpoint  and  techniques  are  presented  to  minimize 
the  number  of  design  variables. 

Although  the  finite  element  method  is  the  most 
extensively  used  discretization  technique,  there  are  some 
reports  of  the  use  of  different  techniques.  Soares  and  Choi 
[13]  presented  a general  formulation  using  the  boundary 
element  method  to  maximize  the  torsional  rigidity  of  a shaft 
or  to  minimize  the  compliance  of  a structure  subject  to  an 
area  constraint. 
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Other  Approaches 

H.  Nooshin  [14]  described  an  algebraic  representation  and 
processing  of  structural  configuration.  A mathematical  system 
introduced  in  the  sequel  consists  of  three  abstract  objects 
referred  to  as  signets,  monads  and  formices,  together  with 
certain  definitions  that  give  rise  to  a number  of  rules  for 
manipulation  of  these  objects.  The  primary  aim  is  to  provide 
a means  for  the  solution  of  problems  related  to  the 
arrangement  and  the  position  of  structural  entities  such  as 
joints,  elements,  loads  and  constraints. 

When  the  problem  size  becomes  very  large  or  the 
constraints  are  implicit  functions  of  the  design  variables, 
applications  of  approximation  concepts  in  structural 
optimization  have  been  utilized.  Kirsch  [15]  presented 
explicit  approximation  models  of  the  optimum  design  problem. 
Ding  and  Gallagher  [16]  presented  a rational  reduced  basis 
reanalysis  and  an  approximation  reanalysis  method  that  is 
based  on  the  size  variable-independence  of  the  equilibrium 
equations  in  matrix  form.  These  techniques  can  be  advantageous 
in  the  structural  optimization  process. 

Mohandas,  Phelps,  and  Ragsdell  [17]  using  fuzzy  goal 
programming  approach  for  the  structural  optimization.  The 
approach  is  based  on  the  assumption  that  the  designer  is  able 
to  specify  target  values  for  various  design  goals.  Instead  of 
formulating  a design  problem  in  terms  of  scalar  optimization 
problem,  a development  of  fuzzy  goal  programming  is  proposed. 
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The  problem  is  solved  by  constructing  a decision  function  as 
a fuzzy  set  using  fuzzy  arithmetic. 

Velazquez  and  Seireg  [18]  presented  a rule-based 
algorithm  for  optimizing  the  shapes.  The  generation  of  the 
optimal  shape  is  governed  by  the  selection  of  an  adequate 
value  of  a critical  stress  for  each  iteration.  This  value 
defines  the  minimum  stress  an  element  should  carry  to  remain 
as  part  of  the  structure . The  new  shape  was  based  on  the 
elimination  and  the  reduction  of  the  elements  at  each  design 
iteration.  It  is  interesting  that  the  rule-based  algorithm 
may  be  converted  to  the  knowledge-based  system  after 
developing  the  rules  needed  for  shape  optimization.  Atrek 
[19]  developed  the  computer  program  SHAPE  to  utilize  the 
elimination  concept  for  shape  optimization  of  continuum 
structures . 

Shape  Optimization  for  Different  Loading  Conditions 
Examples  of  form  synthesis  in  mechanical  structures 
generally  fall  into  one  of  the  following  categories; 

(1)  optimization  for  static  loads  with  constraints 

(2)  optimization  for  dynamic  loads  with  constraints 

(3)  optimization  for  thermoelastic  requirements 

(4)  optimization  for  aeroelastic  requirements 

This  study  deals  with  the  first  two  classes  of  the  problem. 
The  static  optimization  problem  can  be  formulated  as  a linear 
or  nonlinear  programming  problem,  depending  on  the  nature  of 
the  objective  function  and  constraints.  In  such  a problem  the 
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objective  function  is  usually  a differential  function  of  the 
decision  variables.  In  contrast,  in  dynamic  optimization  the 
objective  function  is  usually  a functional  and  frequently  a 
time  integral  from  initial  state  to  the  state  at  another  time. 
The  Static  Problem 

There  are  numerous  published  studies  of  form  synthesis  in 
structural  optimization.  Seireg  [1]  presented  a review  of 
some  illustrative  examples  of  the  use  of  optimization 
techniques  in  the  design  of  mechanical  elements  and  systems. 
Belegundu  and  Arora  [2]  presented  a comprehensive  review  of 
various  mathematical  programming  techniques  used  for 
structural  optimization.  Mathematical  programming  techniques 
[3-7]  were  also  used  in  shape  optimization. 

When  computer-aided  geometric  design  tools  became  more 
and  more  accessible  to  design  engineers,  it  was  a natural 
development  to  couple  CAD  systems  with  structural  optimization 
and  minimum  weight  design  problems  [8-10] . Botkin  and  Bennett 
[11]  used  automatic  mesh  generation,  which  was  based  only  on 
boundary  points  of  the  boundary  design  elements  coupled  with 
adaptive  refinement,  in  the  shape  optimization  of  plate 
structures.  A study  by  Imam  [12]  deals  directly  with  three- 
dimensional  shape  optimization  of  structural  components  by 
using  a concept  of  shape  algorithm. 

Nooshin  [14]  described  an  algebraic  representation  and 
processing  of  structural  configuration.  When  the  problem  size 
is  very  large  or  the  constraints  are  implicit  functions  of  the 
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design  variables,  applications  of  approximation  concepts  in 
structural  optimization  became  necessary.  Kirsch  [15] 
presented  explicit  approximation  models  for  the  optimum  design 
problem.  Ding  and  Gallagher  [16]  presented  a rational  reduced 
basis  reanalysis  and  an  approximation  reanalysis  method  that 
is  based  on  the  size  variable-independence  of  the  equilibrixom 
equations  in  matrix  form. 

There  has  been  increasing  research  activity  in  the  area 
of  design  sensitivity  analysis  in  structural  optimization  [20- 
25].  In  design  sensitivity  analysis,  quantitative  information 
on  how  the  response  of  a structure  is  affected  by  the  change 
in  the  variables  that  define  its  shape  is  obtained  for 
derivatives  of  loads,  stress  and  stiffness  matrices  of  the 
finite  elements  of  the  model.  There  are  different  approaches 
to  design  sensitivity  analysis,  and  any  one  of  them  may  be 
incorporated  into  an  optimality  criterion  or  a mathematical 
programming  method  for  structural  optimization. 

A rule-based  algorithmic  approach  was  applied  for  shape 
optimization  [18,25].  The  generation  of  the  optimal  shape  is 
governed  by  the  selection  of  an  adequate  value  of  a critical 
stress  for  each  iteration.  This  value  defines  the  minimum 
stress  an  element  should  carry  to  remain  as  part  of  the 
structure.  A new  shape  is  determined  from  the  elimination  and 
the  reduction  of  the  elements  at  each  design  iteration.  The 
developed  rule-based  algorithm  may  be  converted  to  a 
knowledge-based  system  after  developing  the  rules  needed  for 
shape  optimization. 
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Form  Synthesis  with  Body  Forces 

Seireg  [26]  used  the  gradient  search  method  to  optimize 
the  profile  of  rotating  disk.  Bhavikatt  and  Ramakrishnan  [27] 
determined  the  optimum  shape  of  rotating  disks  using  a 
nonlinear  programming  approach.  The  profile  is  defined  by  a 
fifth  degree  polynomial  that  is  completely  evaluated  from  the 
boundary  conditions  and  design  variables.  Cheu  [28]  developed 
a master  element  method  for  sensitivity  analysis  in  shape 
optimization  of  axisymmetric  structures.  Park,  Yoo  and  Kwon 
[29]  developed  a design  sensitivity  analysis  for  shape 
optimization  of  an  axisymmetric  turbine  disk  using  the 
boundary  element  method.  De  Silva  [30]  used  a nonlinear 
programming  method  to  minimize  the  weight  of  a circular  disk 
subject  to  specified  behavior  and  side  constraints. 

The  Dynamic  Problem 

Pierson  [31]  gave  a review  of  the  research  on  optimal 
structural  design  under  dynamic  constraints  in  the  early 
1970s.  Frequency  constraints  were  among  the  first  dynamic 
constraints  that  were  considered.  Under  dynamic  conditions  a 
structural  design  problem  generally  consists  of  minimizing  the 
weight  of  the  structure  while  keeping  its  fundamental  natural 
frequency  greater  than  or  equal  to  a specified  value  [32-34]. 
The  dual  problem  is  to  maximize  the  natural  frequency  with 
constrained  weight  [35].  Khot  [36]  investigated  optimization 
of  structures  with  multiple  frequency  constraints,  including 
the  frequencies  being  equal  to  given  values  or  separated  by  a 
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specified  interval  with  a preselected  vibration  mode 
associated  with  the  fundamental  frequency.  The  combination  of 
statical  and  dynamical  constraints  was  also  studied  by  [37-39] 
for  optimizing  the  geometrical  configuration  of  a structure 
subject  to  stress,  displacement  and  frequency  requirements. 

The  frequency  constraints  of  the  structures  influence  its 
vibrational  characteristics  and  consequently  affect  its 
dynamic  response.  The  differential  equations  of  motion  play 
a more  central  role.  The  dynamic  response  constraints  are 
important  in  dealing  with  the  dynamic  loading  optimization. 
In  many  practical  mechanical  and  structural  systems,  the 
transient  dynamic  response  is  of  interest.  The  field  of 
optimization  of  the  dynamic  response  of  mechanical  systems 
[40-44]  generally  deals  with  the  design  of  systems  with 
minimized  velocities,  stresses  or  displacements  during  the 
loading  period. 

The  Objective  of  Study 

The  knowledge-based  shape  optimization  system  is  a target 
that  deserves  further  investigation.  This  study  is  undertaken 
to  develop  a generalized  procedure  for  the  design  of  two- 
dimensional  structures  that  incorporates  elimination  rules, 
mesh  generation,  shape  update,  and  boundary  smoothing.  It  is 
not  an  easy  task  to  incorporate  the  rules  developed  in  this 
study  into  commercial  finite  element  programs.  The  finite 
element  program  used  in  this  study  is  developed  by  the  author 
to  be  suitable  for  the  elimination  concept  of  optimization. 
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Some  special  features  of  the  finite  element  method  are 
included  in  the  program. 


CHAPTER  2 

AUTOMATIC  MESH  GENERATION  AND  NODE  RENUMBERING 

Introduction 

The  finite  element  method  is  a very  powerful  tool  for 
solving  engineering  problems.  For  a complex  system,  it  is 
time-consuming  work  to  prepare  data  for  the  finite  element 
model,  and  mistakes  can  easily  occur  during  the  preparation 
process.  In  order  to  save  time  and  avoid  errors,  an  automatic 
mesh  generator  is  generally  needed.  There  are  numerous 
studies  dealing  with  mesh  generation.  The  transfinite  mapping 
method  [45]  is  commonly  used  in  the  finite  element  generation 
of  the  triangular  or  quadrilateral  boundaries.  The 
flexibility  of  the  transfinite  mapping  method  is  limited  in 
the  mesh  density  adjustment,  and  it  requires  subdivision  of 
the  domain  into  triangular  or  quadrilateral  subdomains  before 
the  mesh  generation. 

From  the  mesh  generation  point  of  view,  a triangular  mesh 
is  easier  to  generate  than  a quadrilateral  mesh.  There  are  a 
wide  variety  of  triangular  mesh  generation  methods  available 
[46,47].  Some  of  these  methods  cannot  be  implemented  within 
the  current  geometric  modeling  theories  and  consequently  do 
not  fit  the  fully  automatic  mesh  generation  requirement. 
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The  recursive  subdivision  method  [48]  is  based  on 
recursively  dividing  the  object  to  be  meshed  into  subdivisions 
and  generating  nodes  along  the  boundaries.  The  method  can  be 
used  to  generate  triangular  or  quadrilateral  elements  within 
any  general  planer  n-sided,  simple-connected  region.  The  mesh 
density  along  the  region  boundaries  can  be  specified  by  the 
user.  Multiple-connected  regions  must  be  reduced  to  simple- 
connected  region  by  connecting  the  boundaries  of  holes  to  the 
outer  boundaries.  Splitting  the  polygon  along  the  "best 
splitting  line"  and  generating  nodes  along  this  new  line, 
based  on  the  mesh  density,  produces  two  subpolygons.  This 
process  is  continued  until  all  subpolygons  have  been  reduced 
to  a simple  polygon  and  these  are  the  finite  elements.  By 
combining  the  recursive  subdivision  method  with  geometric 
modeling  theory,  the  mesh  generator  developed  in  this  study  is 
fully  automatic  with  minimum  interference  by  the  user. 

In  the  procedures  of  shape  optimization,  the  meshing 
routine  should  be  automatic.  The  shape  is  updated  after  each 
design  iteration.  The  user  is  not  expected  to  create  meshes 
at  every  iteration  in  the  optimization  process.  Combinations 
of  lines,  arcs  and  Bezier  splines  are  used  to  smooth  the 
outlines  of  the  shapes.  The  mesh  generation  program  is  needed 
in  every  iteration  of  the  design  optimization  process.  There 
are  some  requirements  that  should  be  satisfied  in  order  to 
achieve  the  design  goal,  and  a fully  automatic  mesh  generator 
should  have  the  following  characteristics: 
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(1)  Data  input  is  easy  and  there  is  a minimum  of  user 
input.  If  the  mesh  generation  needs  lots  of  data  input 
and  preparation,  it  cannot  be  called  fully  automatic, 
because  the  user  must  prepare  too  many  things  before 
starting  the  mesh  generator,  such  as  subdivision  of 
domain.  In  the  shape  optimization  procedure,  the 
boundary  of  structure  is  changed  in  every  iteration.  If 
the  mesh  generator  needs  lots  of  user  input,  the 
optimization  process  has  to  be  interrupted  after  every 
design  iteration  to  prepare  data  for  the  next  iteration. 
If  the  mesh  generator  needs  only  minimum  data  to  generate 
meshes , it  is  simple  and  easy  to  prepare  data  or  to 
access  data.  This  is  also  helpful  for  two  consecutive 
shape  optimization  iterations,  because  it  needs  less 
programming  effort  to  access  and  change  the  data  of 
structure  boundaries . 

(2)  The  density  of  mesh  elements  can  be  controlled  over 
the  whole  region.  The  mesh  size  will  be  decided 
according  to  the  finite  element  calculation  and  the 
design  goal.  This  should  be  set  by  the  user  or  by  the 
program  rules. 

( 3 ) The  mesh  generation  can  handle  multiple  connected 
loops.  The  rule-based  elimination  shape  synthesis  will 
always  create  a hole  inside  the  structure,  and  the 
requirements  for  the  design  also  may  produce  multiple 
holes.  The  mesh  generator  must  have  this  capability  for 
general  two-dimensional  cases. 
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(4)  The  mesh  generator  yields  an  ordering  sequence 
system  for  the  finite  element  processing,  and  the 
characteristics  of  the  mesh  is  suitable  for  the  finite 
element  calculation.  This  is  the  basic  requirement  for 
the  mesh  generation  program. 

In  this  thesis  the  mesh  generation  program  was  developed 
to  fulfill  the  above  requirements  and  was  combined  with  the 
other  modules  for  optimum  shape  design.  The  mesh  generation 
program  was  also  modulized  for  independent  use  only , if 
necessary,  so  it  can  be  used  by  anyone  who  wants  to  generate 
meshes  for  a two-dimensional  multiple  connected  shape. 

Description  of  Structure 

To  represent  the  shape  of  a structure  in  the  computer, 
some  type  of  data  structure  must  be  defined  to  process  the 
boundary  of  shape,  as  shown  in  figure  2—1.  The  data  structure 
must  be  complete  enough  for  most  of  the  conditions  that  may 
occur  in  the  mesh  generation  procedures.  It  should  also  be 
flexible  and  easy  to  process.  There  are  three  basic  geometric 
elements  that  can  be  used  to  describe  the  boundary  of  a shape: 
line,  arc,  and  Bezier  spline.  The  circle  is  a special  case  of 
the  arc.  The  Bezier  spline  is  used  to  smooth  out  the  zigzag 
outline  that  generally  results  in  each  design  iteration  after 
some  elements  are  eliminated  from  the  structure.  If  the  size 
ratio  of  two  connected  elements  is  too  large,  the 
characteristics  of  element  will  be  undesirable  for  finite 
element  calculations.  It  may  cause  a large  distortion  of  the 
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element.  Element  nodes  inside  each  geometric  element  are 
therefore  linearly  placed  between  two  end  points  of  the 
element.  Each  geometric  element  describing  the  shape  of  the 
structure  includes  the  geometric  data  and  mesh  size  data  at 
two  end  points.  The  data  structure  for  every  geometric 
element  of  shape  has  the  following  fields. 

(1)  Element  type;  This  defines  the  type  of  element, 
line,  arc,  or  Bezier  spline. 

(2)  Direction  flag:  For  the  arc  element  this  flag  is 
used  to  indicate  the  direction  of  arc,  counterclockwise 
or  clockwise.  For  the  Bezier  spline  element  this  flag  is 
used  as  the  counter  of  the  number  of  points  of  the 
spline. 

(3)  First  data:  For  the  line  element  this  field  stores 
the  X coordinate  of  the  first  end  point.  It  stores  the 
X coordinate  of  the  center  for  the  arc  element. 

(4)  Second  data:  For  the  line  element  this  field  stores 
the  y coordinate  of  the  first  end  point.  It  stores  the 
y coordinate  of  the  center  for  the  arc  element. 

(5)  Third  data;  For  the  line  element  this  field  stores 
the  X coordinate  of  the  second  end  point.  It  stores  the 
radius  of  the  arc  element. 

(6)  Fourth  data:  For  the  line  element  this  field  stores 
the  y coordinate  of  the  second  end  point.  It  stores  the 
angle  of  the  starting  point  of  the  arc  element. 

(7)  Fifth  data:  It  stores  the  angle  of  the  ending  point 


of  the  arc  element. 
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(8)  Mesh  size:  The  mesh  size  at  the  starting  point  of 

the  geometric  element. 

(9)  Mesh  size:  The  mesh  size  at  the  end  point  of  the 

geometric  element. 

(10)  Array  pointer:  The  pointer  of  array  is  used  to 

store  the  x coordinates  of  all  Bezier  spline  points. 

(11)  Array  pointer:  The  pointer  of  array  is  used  to 

store  the  y coordinates  of  all  Bezier  spline  points. 

(12)  Link  pointer:  Previous  link  list  pointer  of  the 

shape  loop  elements. 

(13)  Link  pointer:  Next  link  list  pointer  of  the  shape 

loop  points. 

The  outline  of  shape  generally  is  a multiple  connected 
loop  and  each  loop  is  created  by  joining  the  geometric 
elements  in  a counterclockwise  sequence.  In  the  shape 

updating  procedure  the  loops  of  shape  outline  are  changed  by 
a Boolean  operation  between  the  outlines  of  shape  and  the 
eliminated  groups.  Every  loop  uses  a double  link  list  data 
structure  to  ensure  the  flexibility  of  the  algorithm. 

For  the  multiple  connected  loops  there  are  some  split 
lines  that  connect  different  loops  into  one  closed  loop  as 
shown  in  figure  2-2.  The  split  line  connection  approach  will 
be  discussed  later  in  this  chapter.  The  mesh  size  of  each 
geometric  element  on  the  boundary  is  specified  at  the  joint 
point  of  two  connected  geometric  elements,  and  the  mesh  size 
has  the  same  value  on  both  jointed  elements  to  ensure  the 
consistency  of  the  element  connection  in  the  loop. 
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Mesh  Node  Generation  of  Shape  Boundary 
A good  mesh  should  have  the  proper  internal  angle  for 
each  corner  of  the  mesh  with  no  distortion  around  the  mesh. 
The  placement  of  mesh  nodes  along  the  geometric  elements  on 
the  boundary  and  split  lines  is  an  important  factor  for  the 
goodness  of  final  mesh.  Proper  spacing  of  nodes  is  a very 
important  aspect  of  the  mesh  generator.  The  nodes  are 

linearly  placed  on  the  element  according  to  the  mesh  size 
specification  of  the  two  end  points  of  the  element.  In  most 
cases  the  element  length  may  not  exactly  match  the  theoretical 
length  needed  for  odd  or  even  subdivisions  of  the  element. 
The  error  between  the  actual  length  and  the  theoretical  length 
should  be  minimized  to  place  the  nodes  inside  the  element. 
First  we  have  to  assume  that  there  are  n internal  nodes  (n+1 
segments)  inside  the  element.  Let  the  number  of  segments  N be 
egual  to  n+1.  The  size  of  each  internal  segment  varies 

linearly  from  the  starting  point  with  mesh  size  to  the  end 
point  with  mesh  size  of  as  shown  in  figure  2-3.  Then  the 
mesh  size  between  two  internal  nodes  is  calculated  by  the 
linear  equation 

Li=m*Ni+b  (2.1) 

where  L.  = mesh  size  at  i-th  subdivision 
= segment  number 

The  parameters  m and  b can  be  solved  from  the  conditions  at 
the  end  points. 

Lg  = m * 1 + b at  starting  point 
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L|^  = in*N  + b at  end  point 
Then  the  parameters  m and  b are 

m = ^ ^b~La  (2.2) 

N-1  n 


b = L - = L - m (2.3) 

an® 

The  length  of  i-th  segment  can  be  expressed  as: 

L.  = i + L, - = m i + (L  -m)  (2.4) 

^ n ® n 

The  summation  of  all  segments  is  equal  to  the  ideal  length  of 
element. 


N 


E i-i 


N 


= E 

i-1 


(2.5) 


This  will  yield  the  ideal  length 


L 


ideal 


N(L^-^Lt,) 

2 


(2.6) 


The  calculated  ideal  length  will  not  be  exactly  equal  to  the 
real  length.  The  ratio  between  the  true  length  and  the  ideal 
length  is  the  scaling  factor  used  to  calculate  real  segment 
size  at  the  internal  nodes. 

Ratio  = (2.7) 

^ideal 


To  calculate  the  i-th  segment  length,  we  have  to  calculate  the 
ideal  length  at  i-th  segment  using  equation  (2.4).  Then  the 
length  at  i-th  segment  is  multiplied  by  the  scaling  factor  in 
equation  (2.7). 
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(Li)  real  = (Li)  ideal  * ^atio  (2.8) 

Using  real  length  for  each  segment,  we  can  place  the  mesh 
nodes  inside  the  geometric  element  on  the  boundaries.  The 
nodes  on  the  boundary  are  also  generated  in  a counterclockwise 
sequence  consistent  with  the  geometric  elements  of  the  shape. 

Visible  Nodes  of  Concave  Loop 
Using  the  split  line  concept  in  the  mesh  generation  the 
visible  point  must  be  decided  before  finding  an  optimum  split 
line.  In  order  to  find  out  all  visible  points  from  a given 
reference  point  Pc  on  the  loop  we  can  define  a visible  cone. 
The  visible  cone  is  defined  by  two  reference  lines.  The  first 
reference  line  uses  the  given  point  Pc  as  the  start  point  and 
the  next  point  to  Pc,  point  Pn,  as  the  end  point.  The  second 
reference  line  also  uses  Pc  as  the  start  point  and  the 
previous  point  to  Pc,  point  Pp,  as  the  end  point. 

The  visible  cone  is  the  area  sweeping  from  reference  line 
one  to  reference  line  two  in  the  counterclockwise  sequence,  as 
shown  in  figure  2-4.  Any  point  on  the  loop  is  not  a visible 
point  when  it  is  located  outside  the  visible  cone.  The  point 
inside  the  visible  cone  is  the  candidate  visible  point  and 
some  further  examination  should  be  done  to  decide  if  it  is 
visible  or  not.  Accordingly,  the  requirements  that  must  be 
satisfied  for  a candidate  visible  point  are: 

(1)  The  candidate  point  must  be  located  between  next 
point  Pn  and  previous  point  Pp  in  the  counterclockwise 
linked  sequence. 
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(2)  The  candidate  point  lies  inside  the  visible  cone. 

( 3 ) The  line  connecting  the  candidate  point  and  the 
reference  point  will  not  cross  any  line  segment  of  the 
loop. 

(4)  There  is  no  point  on  the  loop  too  close  to  the  line 
connecting  the  candidate  point  and  the  reference  point. 
If  there  is  at  least  one  point  that  is  too  close  to  the 
line,  we  should  remove  this  point  from  the  candidate 
list.  Although  it  is  visible  from  the  reference  point, 
a point  too  close  to  the  line  will  cause  a very  small 
dimension  on  a few  meshes.  In  this  case  these  few 
generated  meshes  will  not  fit  the  requirement  of  a well- 
conditioned  mesh.  Therefore,  when  the  program  checks  the 
possible  visible  point,  a proper  tolerance  should  be 
incorporated  in  the  line-crossing  check. 

For  a single  loop  shape,  the  program  checks  all  four 
requirements  and  lists  the  visible  points  from  the  reference 
point  Pc.  The  visible  points  found  here  are  then  used  to 
determine  the  optimum  split  line  for  the  mesh  generation. 
This  will  be  discussed  later  in  the  chapter.  A visible  line 
is  only  the  first  requirement  for  the  best  split  line  to 
subdivision  shape. 

Best  Split  Line  for  Mesh  Generation 
For  one  specific  reference  point  Pc,  we  can  choose  one 
best  visible  line  as  the  best  split  for  this  point.  Every 
point  on  the  loop  has  its  own  best  split  line,  but  only  one  of 
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them  will  be  chosen  as  the  best  split  line  for  the  loop.  In 
this  section  we  will  discuss  the  criteria  for  selecting  the 
best  split  line  for  a point  and  for  a loop.  The  lines 
connecting  the  reference  point  and  its  visible  points  are  the 
candidate  split  lines  from  which  the  best  one  is  to  be 
selected. 

There  are  five  basic  geometrical  properties  for  any  split 
line:  four  internal  angles  and  the  length  of  the  split  line, 
as  shown  in  figure  2-5.  The  value  of  any  of  these  properties 
will  affect  the  mesh  size.  It  is  therefore  necessary  to 
d0termine  what  are  the  best  values  for  both  the  internal 
angles  and  the  length  of  the  split  line.  Another  question 
that  needs  to  be  answered  is  what  kind  of  properties  of  a mesh 
are  needed  in  the  mesh  generation.  For  a triangular  mesh,  the 
equilateral  triangle  is  the  desired  shape,  but  it  is  not  easy 
to  get  equilateral  triangular  mesh  all  over  the  region.  The 
second  choice  is  that  the  largest  internal  angle  of  the 
triangular  mesh  should  be  as  close  as  possible  to  the  internal 
angle  of  the  equilateral  triangle.  This  will  make  all  three 
internal  angles  of  the  triangular  mesh  approach  sixty  degrees. 
Any  one  of  the  internal  angles  of  the  split  line  may  become 
the  internal  angle  of  mesh.  For  a quadrilateral  mesh  no 
internal  angle  should  be  over  one  hundred  eighty  degrees  to 
avoid  distortion.  A rectangular  geometry  with  an  appropriate 
aspect  ratio  is  the  preferred  shape  in  mesh  generation. 
Accordingly,  the  value  of  the  internal  angles  of  a good  split 
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line  should  be  as  close  as  possible  to  achieving  this 
condition.  An  internal  angle  close  to  zero  degrees  or  one 
hundred  eighty  degrees  will  possibly  create  an  ill-conditioned 
mesh,  and  should  be  avoided  in  the  mesh  generation  process. 
The  following  are  several  factors  that  are  included  in  the 
criteria  of  the  split  line  for  a reference  point. 

(1)  The  minimum  of  all  internal  angles  of  the  split  line 
should  be  as  large  as  possible.  This  will  give  the  mesh 
generator  more  flexibility  for  the  next  subdivision 
process.  A small  internal  angle  of  the  split  line  may 
result  in  a mesh  with  small  internal  angle  or  the  split 
line  for  the  next  subdivision  process  will  have  an  even 
smaller  internal  angle. 

(2)  The  length  of  the  split  line  should  not  be  large  in 
comparison  with  the  diagonal  of  the  bounding  box  of  the 
shape.  Proper  tolerance  is  involved  in  the  line  crossing 
test  for  the  visible  point  as  mentioned  before  (see 
figure  2-6).  This  check  reduces  the  importance  of  the 
length  in  the  criteria.  If  there  is  only  computational 
tolerance  when  the  program  tests  the  line  crossing  of  the 
visible  point,  that  point  may  be  very  close  to  the 
visible  line  without  crossing  it.  This  situation  may 
produce  an  unsatisfactory  mesh  with  a very  small  sides. 
Using  the  proper  tolerance,  the  minimum  size  of  any  side 
of  the  mesh  may  be  secured,  guaranteeing  that  the 
generated  mesh  is  acceptable. 
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(3)  The  internal  angles  of  the  split  line  should  keep 
away  from  zero  degrees  or  one  hundred  eighty  degrees. 
The  desired  shape  for  a quadrilateral  meshes  is  a 
rectangle.  The  deviation  from  ninety  degrees  of  every 
internal  angle  of  the  split  line  should  be  minimized. 
Small  deviations  can  keep  the  subdivision  process  on  the 
safe  side. 

From  the  above  three  conditions,  a good  split  line  should 
have  a larger  minimum  of  internal  angles,  a short  length  of 
the  candidate  split  line  for  the  reference  point  and  minimum 
deviation  from  ninety  degrees  for  all  internal  angles.  The 
criteria  for  a good  split  line  can  therefore  be  defined  as 
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where  Cl  : coefficient  of  minimum  internal  angle 
T . : minimum  of  internal  angles 

C2  : coefficient  of  length 
Lg  : length  of  split  line 
C3  : coefficient  of  deviation 
Diag  : diagonal  of  bounding  box  of  shape 
Dev  : deviation  of  internal  angle  from  ninety  degree 
From  equation  (2.9)  we  can  determine  the  best  split  line 
for  a specific  reference  point.  In  order  to  determine  the 
best  split  for  the  loop,  the  global  property  of  the  reference 
point  should  be  included  in  the  consideration.  However,  it  is 
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not  easy  to  define  the  surrounding  circumstances  of  the 
reference  point.  In  the  program  the  number  of  visible  points 
and  the  average  value  of  the  objective  function  of  all  split 
lines  for  the  reference  point  are  used  as  the  global 
properties.  The  more  visible  points  from  the  reference  point, 
the  safer  the  process  is  to  select  it  as  the  best  split  line. 
If  the  average  of  the  objective  function  of  all  split  lines  of 
a reference  point  is  high,  then  the  quality  of  split  line 
using  this  reference  point  is  good.  Consequently  the  best 
split  line  of  this  reference  point  may  be  used  as  the  best 
split  line  of  the  loop.  In  the  program  the  best  split  line  of 
a point  is  selected  first,  then  it  is  compared  with  the  best 
split  line  previously  chosen  in  the  process.  If  they  have  the 
same  objective  function  value,  the  point  with  the  higher 
average  object  function  value  of  all  its  split  lines  is 
selected.  The  procedure  is  repeated  until  all  points  on  the 
loop  are  compared.  After  the  best  split  line  is  picked,  the 
program  subdivides  the  double  link  list  loops  into  two  double 
link  list  loops  by  using  the  split  line  as  the  separator. 

The  triangular  mesh  and  the  quadrilateral  mesh  are  used 
for  the  mesh  generator.  In  the  split  line  algorithm  the 
subdivision  process  is  repeated  until  the  condition  to 
terminate  the  process  is  reached  which  is  the  number  of  nodes 
of  the  loop  are  less  than  or  equal  to  four.  This  allows  the 
loop  to  become  a triangular  or  a quadrilateral  mesh.  In  the 
finite  element  calculation,  the  element  nodes  should  arranged 
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in  same  sequence.  Since  the  nodes  on  the  loop  are  always  kept 
in  a counterclockwise  sequence,  the  final  subdivision  loop 
with  three  and  four  nodes  is  easily  converted  to  an  element  of 
the  mesh.  Initially  when  a loop  is  split  into  two  loops  by  a 
split  line,  the  remaining  loops  will  still  have  more  than  four 
nodes.  Each  of  them  needs  a further  subdivision  process  and 
the  loops  are  therefore  arranged  in  a sequence  like  stack  in 
the  data  structure,  first  in-last  out.  In  this  way  every  loop 
is  sequentially  subdivided  into  meshes. 

If  shape  has  more  than  one  loop,  the  split  line  is  used 
as  a connecting  line  to  join  the  inside  holes  and  the 
outermost  loop.  After  connecting  the  inside  holes  to  the 
outermost  loop,  the  only  loop  remaining  is  the  outermost  loop. 
This  becomes  the  initial  status  for  the  mesh  generation  to 
start  the  subdivision  process  for  the  multiple  loop  connected 
shape.  The  subdivision  procedure  is  then  repeated  until  all 
the  subdivision  loops  have  less  than  four  nodes  and  are 
converted  to  triangular  elements  or  quadrilateral  elements 
depending  on  the  number  of  nodes  of  subdivision.  The  result 
shown  in  figure  2-7  is  an  example  of  a mesh  which  is  generated 
automatically  without  any  constraints. 

Special  Algorithm  for  the  Symmetric  Case 

The  split  line  algorithm  for  mesh  generation  may  create 
a symmetric  mesh  distribution  automatically  for  a symmetric 
shape,  as  shown  before  in  figure  2-7.  This  is,  however,  not 
guaranteed  in  all  situations,  because  it  is  not  easy  to  obtain 
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a unique  value  of  the  objective  function  for  every  node  with 
different  geometric  characteristics.  Even  if  we  can  obtain  a 
unique  objective  function  for  each  node,  there  are  still  some 
cases  where  two  or  more  than  two  nodes  have  exactly  the  same 
geometric  characteristics.  In  such  case  it  is  difficult  to 
decide  which  split  line  should  be  chosen  first  and  therefore 
a problem  is  created  for  the  automatic  symmetric  mesh 
generation.  If  two  split  lines  have  exactly  the  same 
geometric  characteristics  in  the  loop,  we  can  choose  any  one 
of  them  as  the  split  line.  After  that,  the  geometric 
characteristics  of  the  other  one  is  changed.  In  the  symmetric 
case,  there  may  be  two  loops  of  symmetry  about  the  symmetric 
line  with  the  same  geometric  characteristics.  But  inside 
these  two  loops,  there  are  two  split  lines  with  exact  the  same 
characteristics  in  each  loop.  The  program  may  choose  split 
line  Sal  in  loop  one  and  Sb2  in  loop  two  as  shown  in  figure  2- 
8.  The  nonsymmetric  mesh  generation  starts  at  this  point  and 
the  nonsymmetric  characteristics  of  mesh  generation  may  spread 
out  to  the  remaining  regions  that  will  be  subdivided  in  the 
future. 

The  special  algorithm  for  the  symmetric  case  is  developed 
to  generate  a symmetric  mesh  distribution  over  the  entire 
region.  First  the  symmetric  line  is  selected  and  the 
intersections  between  loops  and  the  symmetric  line  are 
calculated.  In  some  cases,  such  as  curved  loops,  the  nodes 
generated  by  segment  generation  along  the  loop  are  not  exactly 
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located  on  the  symmetric  line.  The  number  of  segments  is  kept 
for  the  half  of  loop,  but  the  node  position  is  properly 
adjusted  between  intersections.  The  next  step  is  to  connect 
the  outmost  loop  and  inside  loops  which  have  intersections 
with  symmetric  line.  The  new  shape  to  be  subdivided  becomes 
only  half  of  the  original  shape,  as  shown  in  figure  2-9.  The 
mesh  generator  only  creates  meshes  for  this  half.  After  mesh 
generation  is  finished,  the  program  mirrors  all  meshes  about 
the  line  of  symmetry  and  reverses  the  node  seguence  of  every 
mesh  after  the  mirroring  process  in  order  to  keep  the 
counterclockwise  node  seguence. 

Node  Renumbering 

After  all  meshes  are  generated  by  the  mesh  generator,  the 
node  number  is  created  in  a way  similar  to  the  subdivision 
seguence.  If  we  use  this  node  seguence  for  the  finite  element 
technigues  directly,  it  costs  heavily  in  CPU  time.  The 
bandwidth  and  the  compact  storage  of  the  stiffness  matrix 
would  be  very  large.  A better  way  is  to  renumber  the  node 
sequence  after  mesh  generation  before  passing  the  data  to  the 
finite  element  analysis.  In  this  thesis  the  reversed  Cuthill- 
McKee  ordering  scheme  [49,50]  is  used  to  rearrange  the  node 
sequence  and  reduce  the  bandwidth  of  the  matrix  at  the  same 
time.  A connected  graph  representing  the  node  relationship 
must  be  supplied  for  the  reverse  Cuthill-McKee  scheme.  The 
node  number  associated  with  every  mesh  is  used  to  generate  a 
connected  graph  that  interprets  the  matrix  in  the  finite 
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element  method.  Here  we  will  discuss  the  interface  between 
the  mesh  generator  and  the  node  renumbering  program.  For 
simplicity  we  assume  that  every  node  has  only  one  degree  of 
freedom  such  that  the  node  number  is  the  same  as  the  colvimn 
and  row  positions  in  the  matrix.  The  local  and  global  node 
seguences  for  a guadri lateral  mesh  is  the  basic  information 
needed  to  build  the  graph. 

A graph  G=(V,E)  contains  a set  of  vertices  vl,v2,v3, 
. , together  with  a set  of  edges , where  an  edge  is  a pair 
(vi,vj)  of  vertices  of  V.  A graph  may  be  associated  with  any 
matrix  A.  If  A is  a square  of  order  n,  A.j  are  its  elements. 
The  diagonal  elements  A.,  are  all  different  from  zero. 
Accordingly,  the  graph  is  undirected  and  contains  n labeled 
vertices  vl,v2, . . . ,vn,  and  (Vj,Vj)  is  an  edge  of  the  graph  if 
and  only  if  A.j  not  equal  to  zero.  If  (v-,Vj)  is  an  edge, 
vertices  v.  and  Vj  are  said  to  be  adjacent.  If  Y is  a subset 
of  V,  then  the  adjacent  set  of  Y,  Adj(Y),  is  simply  the  set  of 
nodes  in  G which  are  not  in  Y but  are  adjacent  to  at  least  one 
node  in  Y.  The  edge  (v.  ,Vj)  is  said  to  be  incident  to  the 
vertex  Vj  and  the  vertex  Vj.  The  degree  of  a vertex  Deg(v)  is 
the  number  of  edges  incident  to  it.  From  the  point  of  view  of 
mesh  generation,  the  degree  is  the  number  of  side  of  meshes 
using  the  same  node  v or  number  of  vertices  connected  to  same 
point  V.  The  distance  d(Vj,Vj)  between  two  nodes  v-  and  Vj  in 
the  connected  graph  is  simply  the  length  of  a shortest  path 
jointing  nodes  v-  and  v^.  The  eccentricity  of  a node  v,  e(v) , 
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can  be  defined  as  the  largest  distance  between  v and  any  other 
vertex  of  the  graph.  The  largest  eccentricity  of  any  vertex 
in  a graph  is  defined  as  the  diameter  of  the  graph.  A 
peripheral  vertex  is  one  for  which  the  eccentricity  is  equal 
to  the  diameter  of  the  graph.  A pseudo-peripheral  vertex  v is 
defined  by  the  condition  that  if  u is  any  vertex  for  which 
d(v,u)=e(v),  then  e(v)=e(u).  This  definition  guarantees  that 
the  eccentricity  of  a pseudo-peripheral  vertex  is  large  and 
usually  close  to  the  diameter  of  the  graph. 

The  finite  element  graph  associated  with  a constructed 
system  matrix  A in  which  a row  and  a column  correspond  to  each 
node  in  grid  is  shown  in  figure  2-10.  An  element  of  matrix  Ajj 
is  not  equal  to  zero  if  and  only  if  node  i and  j belong  to  the 
same  finite  element.  There  is  a one-to-one  relationship 
between  vertices  of  the  graph  and  the  nodes  of  the  grid.  Two 
vertices  are  connected  by  an  edge  only  if  the  corresponding 
nodes  belong  to  the  same  finite  element.  In  the  following  we 
will  discuss  how  to  construct  the  finite  element  graph  for 
node  renumbering. 

Assume  that  we  have  two  meshes  with  local  node  sequence 
one  to  four  for  each  mesh,  as  shown  in  figure  2-11.  The  first 
mesh  have  global  node  sequence  twenty,  fifteen,  seven,  and 
thirty-one.  The  element  matrix  is  a full  matrix  since  all 
nodes  are  related.  In  the  global  matrix  only  nodes  of  the 
same  mesh  are  related.  Then  for  each  node  there  is  an 
adjacent  set  associated  with  it  according  to  the  information 
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in  the  global  matrix.  The  first  node  in  the  local  sequence  is 
a global  node  twenty  and  it  is  linked  to  the  node  of  the  local 
sequence  two  to  four,  that  is  global  nodes  fifteen,  twenty, 
and  thirty-one.  The  adjacent  set  of  node  twenty  Adj(V2o) 
becomes  {15,20,31}.  Similarly  we  can  find  an  adjacent  set  for 
node  fifteen,  seven,  and  thirty-one  as  listed  in  table  2-1. 


Table  2-1  The  adjacent  sets  of  the  element  one 


node 

adjacent  set 

20 

V, 

15, 

31 

15 

7, 

20, 

31 

7 

15, 

20, 

31 

31 

7, 

15, 

20 

If  two  meshes  have  a common  side,  connecting  node  fifteen 
and  seven,  the  adjacent  set  of  nodes  of  the  second  mesh  will 
be  created  and  the  adj  acent  set  of  common  node  of  the  two 
meshes  will  be  united  to  update  the  adjacent  set  of  that  node. 
For  example,  the  second  quadrilateral  mesh  has  global  nodes 
fifteen,  six,  twelve,  and  seven.  The  adjacent  sets  of  the 
nodes  are  listed  in  table  2-2. 


Table  2-2  Adjacent  sets  of  element  two 


node 

adjacent  set 

15 

6, 

7, 

12 

6 

7, 

12, 

15 

12 

6, 

7, 

15 

7 

6, 

12, 

15 

These  adjacent  sets  will  be  used  to  construct  the  graph  of  the 
global  matrix. 
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When  two  meshes  are  jointed  by  a side  with  node  fifteen 
and  seven,  the  adjacent  set  of  these  two  nodes  in  the  two 
meshes  should  be  united.  The  result  is  listed  in  table  2-3. 


Table  2-3  Adjacent  sets  of  the  connected  elements 


node 

adjacent  set 

6 

7,  12,  15 

7 

6,  12,  15,  20,  31 

12 

6,  7,  15 

15 

6,  7,  12,  20,  31 

20 

7,  15,  31 

31 

7,  15,  20 

The  process  is  repeated  until  the  adjacent  sets  of  all  nodes 
are  created  and  the  adjacent  sets  of  nodes  shared  by  other 
meshes  are  properly  updated.  At  the  final  step  every  node  has 
its  adjacent  set.  From  the  adjacent  set  we  can  see  the 
bandwidth  of  each  node  in  the  matrix  and  the  largest 
difference  of  sequence  between  the  node  and  its  adjacent. 
This  largest  sequential  difference  may  be  interpreted  as  the 
column  height  or  row  width  of  the  node  in  the  matrix.  If  we 
desire  to  reduce  the  bandwidth,  the  largest  sequential 
distance  between  the  node  and  its  adjacent  should  be 
minimized.  The  adjacent  set  associated  with  every  node 
created  by  mesh  generation  constructs  the  finite  element  graph 
of  the  relationship  between  nodes.  This  graph  can  be  used  by 
reversed  Cuthill-McKee  scheme  to  reduce  bandwidth  of  the 
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The  reversed  CutHill-McKee  ordering  algorithm  for  a 
associated  graph  can  be  described  as  follows. 

(1)  Find  a starting  node  r and  assign  vl  to  r.  The 
effect  of  the  ordering  algorithm  depends  is  highly 
dependent  on  the  choice  of  the  starting  node.  The 
objective  is  to  find  a pair  of  nodes  which  are  at  a 
maximum  or  near  a maximum  distance  apart,  pseudo- 
peripheral nodes.  Here  we  have  to  introduce  the  rooted 
level  structure  of  a node.  A level  structure  Lq,  , . . . , 

with  m+1  levels  is  obtained  when  subsets  are  defined 
in  such  a way  that: 

Adj(L.)  c Lj.,  u L.^i,  0 < i < m 
Adj(Lo)  c L, 

Adj(I^)  c Vi 

If  the  node  r is  chosen  as  root,  m is  the  length  of  the 
level  structure  and  is  egual  to  the  eccentricity  e(r). 
To  find  a pseudo-peripheral  node  is  to  find  a rooted 
level  structure  which  has  largest  length  m or  largest 
e(v)  . 

(2)  For  i=l,...,N,  find  all  the  unrenumbered  neighbors 
of  the  node  vi,  and  number  them  in  increasing  order  of 
degree.  N is  the  total  number  of  nodes  to  be  renumbered. 
The  degree  of  a node  vi  is  the  number  of  edges  incident 
to  it  or  the  number  of  nodes  in  th  adjacent  set. 
Selecting  a node  with  a smaller  degree  as  having  higher 
priority  to  be  numbered  first  in  the  node  renumbering 
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process.  This  will  reduce  the  difference  of  the  node 
number  of  a node  vi  and  its  adjacent  nodes  and  will 
consequently  reduce  the  bandwidth, 

(3)  The  reverse  CutHill-Mckee  ordering  is  given  by  u^, 
U2,....,  u^  where  u.  = for  i=l,...,N.  In  the  study 
of  the  profile  methods  by  George  [50] , he  discovered  that 
the  ordering  obtained  by  reversing  the  CutHill-McKee 
ordering  often  turns  out  to  be  much  superior  to  the 
original  ordering  in  terms  of  profile  reduction,  although 
the  bandwidth  remains  unchanged. 

In  simple  test  case,  such  as  a square  region,  the  reverse 
CutHill-McKee  algorithm  sometimes  does  not  produce  the  optimum 
ordering  sequence,  but  it  is  generally  satisfactory.  In 
complex  mesh  generation  cases  with  multiple  connected  holes, 
the  algorithm  does  generate  a very  good  ordering  sequence. 
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Figure  2-2  Boundary  description  of  multiple 
connected  loops 
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Figure  2-7  Illustration  of  automatic  mesh 
generation 


with  same  geometric  characteristics 
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Figure  2-9  Special  algorithm  for  symmetric  case 
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Figure  2-11  Two  connected  finite  elements 


CHAPTER  3 

FINITE  ELEMENT  METHOD 
Introduction 

The  finite  element  method  [51-54]  is  used  in  this  study 
to  calculate  the  stress  distribution  needed  for  shape 
optimization  of  structures  subjected  to  static  loading, 
inertial  loading,  and  impact.  The  result  of  the  finite 
element  analysis  is  utilized  in  the  optimization  procedure. 
In  order  to  evaluate  the  result  quickly,  a postprocessor  is 
also  developed  to  display  the  stress  contours,  stress 
distribution,  and  displacements.  The  finite  elements  are 
generated  by  the  mesh  generation  program  automatically  at 
every  design  iteration.  The  shape,  node  numbering,  and  mesh 
distribution  can  be  different  for  each  iteration  since  it  is 
difficult  for  the  user  to  update  the  boundary  conditions  and 
loads  every  time  the  design  process  starts.  Some  special 
algorithms  are  also  developed  to  handle  the  boundary 
constraints,  the  point  loads,  and  the  surface  pressure 
automatically.  The  design  iterations  can  therefore  proceed 
from  the  beginning  to  the  end  without  user  interference. 

Outline  of  Finite  Element  Method 
In  finding  an  approximate  solution  of  the  continuum 
problem  by  the  direct  methods  of  variational  calculus,  it  is 
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necessary  to  assume  some  admissible  distribution  of  the 
required  quantities  in  terms  of  the  functions  which  are 
defined  throughout  the  region  of  the  continuum.  This  may 
transform  the  continuum  with  its  infinite  degrees  of  freedom 
into  a substitute  finite  degrees  of  freedom  system  whose 
characteristics  approximate  those  of  the  original  model. 
Generating  suitable  approximating  functions  to  satisfy  the 
prescribed  boundary  conditions  is  simple  when  the  geometry  of 
the  continuum  of  interest  is  simple.  But  the  task  becomes 
progressively  more  difficult  as  the  region  becomes  more 
complex  in  shape.  In  the  finite  element  method,  we  may  assume 
that  the  continuum  is  divided  into  a number  of  simple  shape 
sub-regions  of  elements  of  finite  size.  The  distribution  of 
the  required  quantities  is  assumed  to  be  relatively  simple 
within  these  elements.  By  using  a sufficient  number  of 
elements,  we  may  obtain  an  acceptable  approximate  solution  to 
represent  the  overall  real  situation. 

Within  each  sub-region  or  element,  the  distribution  of 
the  required  functions,  such  as  displacements  and  stresses, 
are  assumed  in  a suitable  form  in  terms  of  some  parameters. 
Generally  the  form  of  the  function  is  a polynomial  and  the 
parameters  are  the  values  of  nodal  coordinates.  On  the  basis 
of  these  assumed  distributions,  the  element  properties,  such 
as  strain  energy,  may  be  computed  and  the  summation  over  all 
elements  yields  the  total  strain  energy  stored  in  the  system. 
By  applying  an  appropriate  variational  principle,  we  may 
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obtain  sets  of  algebraic  equations  in  which  the  nodal 
displacement  values  of  the  required  functions  appear  as  the 
basic  unknowns  to  be  determined. 

The  procedure  of  solving  a stress  analysis  problem  by  the 
finite  element  method  can  be  roughly  divided  into  the 
following  steps; 

(1)  Divide  the  continuum  into  a finite  number  of 
elements  of  simple  geometry.  The  number,  size,  type  and 
arrangement  depend  on  the  problem  to  be  solved.  When 
discretizing  the  continuum,  two  distinct  elements  can 
have  common  points  on  the  common  boundary  and  no 
overlapping  is  allowed.  Besides  the  assembled  elements 
should  approximate  the  geometry  of  the  real  domain  as 
closely  as  possible.  This  is  done  by  the  mesh  generator. 

(2)  Select  the  displacement  model.  In  general, 
polynomial  interpolation  functions  are  used  in  the 
solution  model.  The  number  of  terms  used  depends  on  the 
type  of  the  problem  and  the  precision  requirement.  The 
approximate  function  and  certain  of  its  derivatives  must 
be  continuous  within  the  element  and  there  must  be 
compatibility  between  adjacent  elements. 

(3)  Formulate  the  system  equilibrium  equations.  Using 
the  interpolation  functions  assumed  above,  the  total 
strain  energy  of  the  discrete  approximation  to  the 
continuum  may  be  determined  from  the  summation  of  the 
strain  distribution  within  individual  elements.  Then  we 
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can  determine  the  equilibrium  equations  for  the  nodes  of 
the  discretized  continuum  in  terms  of  element 
contribution . 

(4)  Solve  the  system  of  equilibrium  equations  for  the 
nodal  displacements.  After  the  incorporation  of  the 
boundary  conditions,  a solution  for  the  unknowns  of  nodal 
displacement  can  be  obtained.  The  numerical  procedures 
used  in  solving  the  simultaneous  equations  may  be  the 
algorithm  which  is  derived  by  direct  factorization  or  by 
iteration.  There  are  many  methods  commonly  used  in 
finite  element  programs,  such  as  band  solver,  envelop 
solver  and  frontal  solver.  Selection  of  the  solver 
depends  on  the  hardware  used  for  the  program  and  on  the 
requirement  of  speed.  For  the  small  scale  problems,  the 
in-core  solution  by  band  solver  or  skyline  solver  is 
efficient.  If  the  memory  size  is  restricted  by  hardware 
or  if  the  scale  of  the  problem  is  large,  the  active  band 
solver  or  frontal  solver  is  a good  selection  by  using  the 
large  space  in  the  secondary  storage  device  to  store  the 
matrix  in  blocks.  In  this  study,  the  in-core  skyline 
method  is  used  as  the  solver  in  the  finite  element 
procedures . 

(5)  Determine  the  element  strains  and  stresses.  Once  the 
nodal  displacements  are  solved,  the  element  strains  are 
readily  calculated  from  the  displacement  based  model 
using  strain-displacement  relations.  The  stresses  are 
then  obtained  by  means  of  Hooke's  law. 
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Select  Elements  and  Displacement  Functions 

In  the  mesh  generation  triangular  and  quadrilateral 
elements  are  used.  The  constant-strain  triangular  and  four- 
node  quadrilateral  elements  are  selected  for  the  finite 
element  calculation.  It  is  easy  to  derive  the  shape  function 
for  triangular  and  rectangular  elements  in  terms  of  a global 
(x,y)  coordinate  system.  But  it  is  not  as  easy  for  a 
quadrilateral  element.  Accordingly  both  elements  use 
isoparameteric  formulations.  A finite  element  is  said  to  be 
isoparameteric  if  the  same  interpolation  formulas  define  both 
the  geometry  and  the  displacement  shape  function.  Such 
elements  can  satisfy  both  the  geometric  and  displacement 
compatibility  conditions.  The  interpolation  function  is 
constructed  such  that  the  coordinates  of  any  node  along  an 
interelement  bovindary  are  the  same  regardless  of  from  which 
element  the  node  is  approached.  In  this  case  we  can  insure 
that  coordinates  are  compatible  along  that  boundary.  Another 
reason  using  isoparametric  formulas  is  that  calculations  of 
element  matrices  need  the  evaluation  of  many  integrals.  It 
may  be  difficult  or  impossible  to  calculate  these  integral 
analytically.  The  aid  of  the  computer  makes  it  possible  that 
the  integration  is  accomplished  numerically. 

The  shape  functions  are  derived  on  the  basis  of 
displacement  relationships.  For  the  plane  stress  case,  the 
generic  displacements  at  any  point  within  the  element  may  be 
expressed  as  the  column  vector  {U} : 
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{U}  = [ U,  V (3.1) 
where  u and  v are  translations  in  the  s and  t directions.  The 
displacements  within  an  element  must  be  uniquely  defined  by 
the  nodal  values  of  the  element.  The  representation  of 
displacement  u may  be  assumed  to  be  a linear  polynomial 

u = a^+a2S+a3t+  ...  = [ 1st  ....  ] (a)  = [P]  {a}  (3.2) 
where  (a)  represents  the  generalized  variables  of  the  element 
and  [P]  is  the  interpolation  function  which  is  generally  a 
polynomial.  The  number  of  terms  in  the  polynomial  is  decided 
by  the  number  of  node  and  degree  of  freedom  in  the  element. 
Assuming  n terms  in  the  polynomial,  the  n constants  can  be 
evaluated  by  solving  n simultaneous  equations  which  will  arise 
if  the  nodal  coordinates  are  inserted  and  the  displacements 
equated  to  the  appropriate  nodal  displacements.  For  example 


Ui  = a, 

+ 

^2 

+ 

^3 

^ • • • • 

= [ 

1 

Si 

• • • • 

] 

(a) 

U2  = a^ 

• 

+ 

^2 

+ 

^3 

^ • • • • 

= [ 

1 

S2 

^^2  • • • • 

] 

(a) 

• 

= ^1 

+ 

^2 

Si 

+ 

^3 

tl 

^ • • • • 

= [ 

1 

Sn 

tn  .... 

] 

(a) 

In  matrix  form,  the  above  equations  may  be  expressed  in  the 
form 

{U,}  = [PJ  (a)  (3.3) 

The  unknowns  in  (a)  are  obtained  by 

{a}  = (3.4) 

In  order  to  define  (a)  uniquely  in  terms  of  {U^} , the  matrix 
[P^]  in  equation  (3.4)  should  not  be  singular.  This  depends  on 
the  choice  of  the  polynomial  basis  and  the  coordinates  of  the 
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nodes  of  the  elements . Substituting  equation  (3.4)  into 
equation  (3.2),  we  may  get 

u = [P]  {a}  = [P]  (3.5) 

or 

u = [N]  {U^>  (3.6) 

where 


[N]  = [P]  [Pj-i  (3.7) 

Similarly  the  displacement  v can  be  expressed  as 

V = [N]  {V^}  (3.8) 
For  the  reference  element  of  a three-node  constant-strain 
triangular  element  as  shown  in  figure  3-1,  the  polynomial  is 
selected  as 

u = a^  + a2S  + a3t  = [ 1st  ] (a)  = [P]  {a}  (3.9) 
By  substituting  the  node  coordinates  into  equation  (3.2),  we 
obtain  [P„]  as 


10  0 
110 
10  1 


(3.10) 


After  inverting  matrix  [P^]  in  equation  (3.10),  we  get 


[Pn] 


10  0 
-110 
-10  1 


(3.11) 


We  may  substitute  the  equations  (3.9)  and  (3.11)  into  equation 
(3.5)  to  obtain  the  shape  functions  for  the  triangular  element 

[N]  = [ N1  N2  N3  ] (3.12) 


52 


where 

Nl(s,t)  = 1-s-t,  N2(s,t)  = s,  N3(s,t)  = t (3.13) 
The  polynomial  for  four-node  isoparameteric  quadrilateral 
element,  as  shown  in  figure  3-2  is  expressed  in  the  form 

[P]  =a^  + a2S  + a3t  + a^st  = [ Istst  ] (a)  (3.14) 
Substituting  the  node  coordinates  of  the  reference  element 
into  equation  (3.14),  we  get 


[Pn] 


'1  -1  -1  1' 
1 1-1-1 
1111 
1-1  1-1 


The  inversed  matrix  of  [P^]  is 


[Pn] 


■ 1 1 1 1' 
1-1  1 1-1 
4-1-1  1 1 

1-1  1-1 


(3.15) 


(3.16) 


In  a similar  way,  we  obtain  the  shape  functions  for  the 
quadrilateral  from 

[ N ] = [ N1  N2  N3  N4  ] (3.17) 

where 


N1  (s,  t)  =4  (1-s)  (1-t)  , N2  (s,  t)  (1+s)  (1-t) 

4 4 (3.18) 

N3  (s,  t)  (1+s)  (1+t)  , N4  (s,  t)  (1-s)  (1+t) 

4 4 

Formulation  of  System  Equilibrium  Equations 
The  solution  of  problems  in  the  theory  of  elasticity  can 
be  obtained  by  one  of  two  methods.  One  can  directly  solve  the 
governing  equations  for  the  specific  boundary  conditions  or 
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one  can  minimize  an  integral  quantity  related  to  the  internal 
and  external  work  done  by  the  stress  components  or  the  applied 
loads.  The  finite  element  formulation  of  elasticity  problems 
utilizes  the  second  approach.  There  are  three  basic 
formulations  in  finite  elements,  displacement  based 
formulation,  stress (or  force)  based  formulation  and  hybrid 
formulation.  If  we  start  by  selecting  the  displacement 
equations  that  satisfy  the  displacement  boundary  conditions, 
we  minimize  the  potential  energy  in  the  system.  The 

displacement  field  must  satisfy  the  continuity  condition.  If 
we  start  by  selecting  stress  relationships  that  satisfy  the 
stress  bo\mdary  conditions,  we  should  minimize  the 

complementary  energy.  The  stress  field  must  satisfy  the 
equilibrium  requirement.  The  displacement  and  stress  fields 
are  utilized  simultaneously  in  the  hybrid  method.  The 
displacement  based  formulation  is  commonly  used  and  is  the 
method  applied  in  this  study.  The  cases  involved  in  this 
study  are  two-dimensional  elasticity  problems.  The  stresses 
and  strains  used  in  the  following  formulation  are  the  only 
quantities  for  the  stress  case. 

The  strain  energy  for  a differential  element  of  volume  dV 
is  given  by 

6n  = {a}  - ^{CgHa}  (3.19) 

where  (e)  is  the  total  strain  and  the  {e^^}  is  an  initial 
strain.  611  is  called  the  strain  energy  density  and  the  total 
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strain  energy  is  obtained  by  integrating  over  the  volume 
giving 

n = |A({6}T{a}-{eo}T{o})dV  (3.20) 

The  column  vectors  {e}  and  (a)  depend  on  the  problem  being 
solved.  For  this  study,  the  vectors  for  the  case  of  plane 
stress  are 

{6}  = [6^  €yy  Yxy]^  (3.21) 

{o}  = [a^  Oyy  (3.22) 

In  the  theory  of  elasticity  there  are  two  important 
relationships:  Hooke's  law,  which  relates  the  stress  and 

strain  components,  and  the  strain-displacement  relationships. 
Hooke's  law  has  the  general  form 

{o}  = [D]  {€}  - [D]  {€o>  (3.23) 

where  the  [D]  matrix  contains  the  elastic  constants  of  the 
material.  For  plane  stress  the  [D]  matrix  for  isotropic 
material  is 


[D] 


1 V 


E 

l-v2 


V 1 
0 0 


0 

0 

1-v 

2 


(3.24) 


where  E is  the  Young's  modulus  and  v is  the  Poisson's  ratio. 
The  strain-displacement  relationships  are 


e 


XX 


dx  ' ^ dy  ' 


du  ^ ^ 
dy  dx 


(3.25) 
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where  the  u and  v are  displacement  components  in  the  x and  y 
coordinate  directions.  These  displacements  are  written  in 
terms  of  the  nodal  values  in  equations  (3.6)  and  (3.8).  The 
application  of  equation  (3.25)  allows  the  strain  vector  {e>  to 
be  written  in  terms  of  the  nodal  displacements  {U>.  The 
general  form  of  this  relationship  is 

{€}  = [B]  {U}  (3.26) 

where  the  [B]  strain  matrix  is  derived  by  performing  the 
differentiation  of  shape  function  [N]  which  will  be  discussed 
later. 

Since  the  finite  element  formulation  used  in  this  study 
depends  on  the  minimization  of  the  potential  energy,  the 
requirement  is  that  the  displacement  equations  selected  must 
satisfy  the  displacement  boundary  conditions.  The  total 
potential  energy  of  an  elastic  system  can  be  separated  into 
two  components,  one  resulting  from  the  strain  energy  in  the 
body  and  the  other  is  related  to  the  potential  energy  of  the 
internal  and  applied  loads.  The  total  potential  energy  IIp  can 
be  written  as 

rip  = n + Wp  (3.27) 
where  II  is  the  strain  energy  and  Wp  is  the  potential  of  the 
applied  loads.  The  work  done  by  the  loads  is  the  negative  of 
their  potential  energy  or 

W = - Wp  (3.28) 
Combining  equations  (3.27)  and  (3.28)  yields 

IIp  = n - W 


(3.29) 
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Since  the  region  of  interest  is  subdivided  into  a number  of 
elements,  the  total  potential  energy  may  be  written  as  a 
summation  of  the  element  quantities. 

lip  = ^ (n«-W®)  (3.30) 

i=l  i=l 

Utilizing  equations  (3.20)  and  (3.26),  the  strain  energy  for 
a single  element  can  be  written  as 

n«  = |^|({U}MB]MD]  [B]  [U]  - 2{U}'^[B]'^[D]  {€„} 

{€o}^[D]  {€o})dv 

The  last  term  in  equation  (3.31)  is  not  a function  of  the 
nodal  value  (U) . Therefore,  it  has  no  influence  on  the 
minimization  process. 

The  work  done  by  the  applied  loads  can  be  separated  into 
three  distinct  parts:  the  work  due  to  the  concentrated  loads 
Wc,  the  work  resulting  from  the  pressure  on  the  outside 
surface,  and  the  work  done  by  the  body  forces.  The  work  done 
by  concentrated  loads  is  simply  the  value  of  the  force 
multiplied  by  the  distance  through  which  it  acts.  Denoting 
the  nodal  forces  by  {P>  and  the  nodal  displacements  by  (U), 
the  work  is  given  by  the  matrix  product 

= {U}T{P}  = {P}T{U}  (3.32) 

This  definition  assume  that  the  forces  have  been  resolved  into 
components  parallel  to  the  displacement  components. 

The  work  done  by  body  forces  B„  and  B„  within  the  element 

A y 

is  given  by 
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Wb  = J ( U + V By  ) dv  (3.33) 

where  u and  v are  the  x and  y components  of  the  displacement 
within  the  element . Combining  equations  (3.1),  (3.6),  (3.8), 
and  (3.33),  the  work  may  be  expressed  in  the  form 

Wb  = f {U}T[N]’’ 

J V 

The  work  done  by  the  distributed  loads  that  act  on  the  surface 
is 


dv  = f {U}'^[N]'^{B^}dv  (3.34) 

J V 


= 


=/ < 


u + V ) ds 


(3.35) 


where  u and  v are  the  displacement  components,  and  P^  and  P„ 
are  the  stress  components  parallel  to  the  x and  y coordinate 
directions . Combining  equations  (3.1),  (3.6),  (3.8),  and 
(3.35),  the  work  Wp  is  rewritten  in  the  form 


Wp  = 


P 

= f {U}’^[N]’’  ""  ds  = f {U}T[N]T{PJ  ds  (3.36) 

Js  [PyJ  Js  ® 


The  combination  of  ( 3 . 30 ) , ( 3 . 31 ) , ( 3 . 32 ) , ( 3 . 34 ) , and  (3.36) 
yields  the  following  equation  for  the  total  potential  energy. 


IIp  = E[/^-|{U}MB]T[D]  [B]  [U]dv' 


-E 

-E 


J^{U}T[B]T[D]  {€o}dv]  - 5i[|^{U}T{B^}dv 

f {U}T[N]  {P3}dsl  - {U}'"{P} 

J s 


(3.37) 


To  minimize  IIp,  we  differentiate  equation  (3.37)  with  respect 


to  (U)  and  set  it  equal  to  zero. 
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^ =0 


6{U} 


= E 
-E 


/^[B]^[D]  [B]dv{U}j  - (e„)dv]  (3.38) 

'/  lN]^{B,}dv|  - ]■  [Nl'f{P3)ds] 


- P 


We  may  write  the  system  equilibrium  equation  in  the  form 

[K]  {U}  = {Pj,}  + {F^}  + {T3}  + {P)  = {R}  (3.39) 

where 


[K]  = 5;^  [k]  [B]'^[D]  [B]dv  (3.40) 


{T,}  =5^{t3}  [N]MP3}ds  (3.41) 


{F,}  = = E /^[N]^{B,}dv  (3.42) 

The  characteristic  of  the  dynamic  or  wave  propagation 
problem  is  that  the  response  of  the  system  under  consideration 
changes  with  time.  Equation  (3.39)  is  a statement  of  the 
static  equilibrium  of  the  element  assemblage.  In  this 
equilibrium  consideration,  the  applied  forces  may  vary  with 
time,  in  which  case  the  displacements  also  vary  with  time. 
Equation  (3.39)  is  a statement  of  equilibrium  for  any  specific 
point  in  time.  For  the  impulsive  loading  case,  the  inertia 
forces  need  to  be  considered.  Then  a truly  dynamic  problem 
needs  to  be  solved.  The  principle  of  D'Alembert  asserts  that 
any  law  of  statics  becomes  transformed  into  a law  of  kinetics 
if  the  driving  forces  are  augmented  by  inertia  forces.  Using 
this  principle,  the  principle  of  virtual  work  may  extend  to 
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kinetics.  We  can  simply  include  the  element  inertia  forces  as 
part  of  the  body  forces  in  equation  (3.39).  Assuming  that  the 
element  velocities  and  accelerations  are  approximated  in  the 
same  way  as  the  element  displacements,  the  contribution  from 
total  body  forces  to  the  load  vector  {F^>  is 


F,  = 5^  {f,}  = ( {B^-p  [N]  {U}  ) dv 


(3.43) 


{F,}  = Y:  liv)  = E 

-^L  [N]  p [N]'^[N]  {U}dv 

The  {U}  is  the  nodal  point  accelerations  and  p is  the  mass 
density  of  element.  Substituting  equation  (3.44)  into 
equation  (3.39)  we  may  obtain  the  equilibrium  equation  for  the 
dynamic  case 


[K] {U}  + [M] {U}  = {R} 
where  the  mass  matrix  for  the  structure, 

[M]  ='£[m]  =Y^  f^p  [N]  T [N]  dv 

System  equilibrium  equation  (3.39)  and  (3.45)  are  used  to 
solve  the  displacements  for  the  static  or  dynamics  cases. 
Evaluation  of  the  Element  Matrices 
In  the  parametric  formulation  the  element  displacements 
are  interpolated  in  the  same  way  as  the  geometry 


(3.45) 


(3.46) 
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U = , V = (3-47) 

i-l  i-l 

where  u and  v are  the  local  displacements  at  any  point  of  the 
element  and  Uj  and  are  the  corresponding  element 

displacements  at  its  nodes.  In  order  to  calculate  the 
stiffness  matrix  of  an  element,  we  need  to  calculate  the 
strain-displacement  transformation  matrix.  The  element 
strains  are  obtained  in  terms  of  element  displacements  with 
respect  to  the  local  coordinates.  Because  the  element 
displacements  are  defined  in  the  natural  coordinate  system 
using  equation  (3.47),  we  need  to  relate  the  x and  y 
derivatives  to  the  s and  t derivatives.  The  derivatives  of 
shape  functions  frequently  appear  in  the  integrands.  If  we 
wish  to  perform  integrations  over  the  undistorted  element,  the 
integrals  must  be  transformed  into  ones  that  contains  only  s 
and  t.  In  an  isoparametric  representation  we  may  use  the 
following  representation  for  the  x and  y coordinates  within  an 
element 

n n 

x = J]NiXi,  y = J^N^y^  (3.48) 

i=l  i=l 

Equations  (3.48)  are  very  useful  in  that  they  allow  us  to  map 
a point  in  the  local  coordinate  system  (s,t)  to  a point  in 
(x,y)  coordinate  system.  Therefore  equations  (3.47)  and 
(3.48)  are  said  to  provide  an  isoparametric  mapping  from  the 
undistorted  element  to  the  distorted  element.  This  will 
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enable  us  to  take  a distorted  element  and  map  it  into  an 
undistorted  element  for  the  purpose  of  evaluating  integral. 
Note  that  equation  (3.48)  if  of  the  form 

X = f^(s,t)  y = fy(s,t) 

where  f^  and  f„  are  the  function  of  mapping.  The  inverse 
relationship  is 

s = f3(x,y)  t = ft(x,y) 

We  need  the  derivatives  d/dx  and  d/dy  and  it  seems  natural  to 
use  the  chain  rule  in  the  following  form 


_d_  ^ 

d ds  ^ d dt 

(3.49a) 

dx 

ds  dx 

dt  dx 

_d_  ^ 

d ds  ^ d dt 

(3.49b) 

dy 

ds  dy 

dt  dy 

To  evaluate  d/dx  and  d/dy  in  equation  (3.49),  we  need  to 
calculate  ds/dx,  dt/dx,  ds/dy,  and  dt/dy  which  implies  that 
the  explicit  inverse  relationships  in  equation  (3.49)  need  to 
be  evaluated.  These  inverse  relationships  are  in  general 
difficult  to  establish  explicitly,  and  it  is  necessary  to 
evaluate  the  required  derivatives  in  the  following  way.  Using 
the  chain  rule,  we  have 


d] 

' dx 

dyl 

[ d] 

ds 

ds 

ds 

dx 

d 

dx 

dy 

d 

dt 

dt 

dt 

. 

The  2x2  matrix  on  the  right-hand  side  is  known  as  the  Jacobian 
matrix  and  is  denoted  by  J,  or 
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[J] 


' dx 

ds 

ds 

dx 

dt 

at 

(3.51) 


The  X and  y coordinates  within  an  element  is  represented 
by  shape  functions  in  equation  (3.48).  We  may  then  evaluate 
Jacobian  matrix  as 


[dx 


at 

Then  the  inversion  of 

[J]  = 


Jacobian  matrix  is 


ds 

dt 

ay 

_ dy] 

dx 

dx 

1 

ds 

dt 

ds 

dt 

detJ 

dx 

dx 

dy 

dy. 

ds 

dt 

(3.52) 


If  equation  (3.50)  is  premultiplied  by  J'\  we  obtain  the 
desired  result 


a 

' dx 

ay] 

-1 

a 

dx 

ds 

as 

as 

a 

dx 

ay 

a 

.dy. 

dt 

at. 

at 

The  strain  at  a point  within  an  element  may  be  calculated  by 
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{€}  = 


5u 

dx 

dv 

dy 

du_^  dv 
dy  dx 


r=< 


E dx"^^ 

i=l 

' dNj^ 
dx 

0 

A dN, 

n 

dN. 

S ay  ''' 

i-l 

0 

dY 

n 

dy  ^dx 

dN^ 

dNj^ 

E 

U=i 

dy 

dx\ 

\u, 


V 


[ (3.55) 


Then  the  strain  matrix  [B]  for  the  triangular  element  is 

[B]  = [ Bi  B2  Bj  ] (3.56) 

Similarly,  the  strain  matrix  for  the  quadrilateral 

[B]  = [ Bi  B2  B3  B^  ] (3.57) 

where 


B. 


dN^ 

dx 


0 


dNj^  dN^ 
dy  dx 


1=1, . . . ,n 


(3.58) 


Using  equation  (3.58),  we  can  transform  the  [B]  matrix 
from  the  (x,y)  coordinate  system  into  the  (s,t)  coordinate 
system.  To  integrate  the  element  stiffness  matrix 
numerically,  the  limit  of  each  integration  must  be  changed  to 
become  -1  and  1.  The  infinitesimal  volume  for  a constant 
thickness  structure  dv  = tdA  = tdxdy  must  be  replaced  by  an 
appropriate  expression  in  terms  of  ds  and  dt.  Vector  r 
locates  a generic  point  in  the  Cartesian  coordinates  x and  y 
as  follow 

r=x+y=xi+yj 

The  rate  of  change  of  r with  respect  to  s is 
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3r  _ ^ dy  . 

ds  ds  ds^ 


(3.59) 


Also  the  rate  of  change  of  r with  respect  to  t is 


o.eo, 

When  equation  (3.59)  is  multiplied  by  ds  and  equation  (3.60) 
is  multiplied  by  dt,  they  form  two  adjacent  sides  of  the 
infinitesimal  parallelogram  of  area  dA.  This  area  may  be 
determined  from  the  following  vector  triple  product 


^ ^ If  • * 


(3.61) 


Substitution  of  Eqs.  (3.59)  and  (3.60)  into  Eg.  (3.61)  gives 


dA  = 


I dx  dy  _ dy  dx. 
\ ds  0t  3s  3t 


|3s3t 


(3.62) 


The  expression  in  the  parentheses  of  Eg.  (3.62)  may  be  written 
as  a determinant.  That  is 


dA  = 


dx  dy 


ds  ds 
dx  dy 


dsdt  = |j|dsdt 


3t  3t 


(3.63) 


in  which  J is  the  Jacobian  matrix  and  |J|  is  its  determinant. 
The  element  stiffness  matrix  may  be  expressed  as 

[k]  = j’^[B]'^[D]  [B]dv  = [B]^[D]  [B]  |j|dsdt  (3.64) 

The  formulation  we  discussed  before  may  be  used  for  another 
approach.  An  axisymmetric  solid  is  defined  as  a three- 
dimensional  body  that  may  be  developed  by  rotating  a planar 
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section  about  an  axis.  This  type  of  body  sometimes  is  called 
a solid  of  revolution.  We  assume  that  a body  is  axisymmetric 
with  respect  to  the  z axis  and  that  a typical  finite  element 
is  a circular  ring.  If  the  loads  on  an  axisymmetric  solid  are 
also  axisymmetric,  we  may  analyze  a representative  cross 
section  as  if  it  were  a two-dimensional  problem.  For  any 
point  on  an  axisymmetrically  loading  ring  element,  the  generic 
elements  are  also  defined  as 

{U}  = [ u,  V 

Translations  of  u and  v occur  in  the  r and  z directions.  With 
axisymmetric  loads,  the  translations  w in  the  0 direction  is 
zero,  and  the  shearing  strains  and  y^g  are  also  zero. 
However,  four  types  of  strain  that  are  not  zero  may  be 
expressed  as 

{e)^  = [ e,  6g  ] (3.65) 
Relationships  between  the  strains  and  generic  displacements  in 
equation  (3.65)  are  seen  to  be 


du 

3r 

dv 


€e  = 


_ 271  (r +u)  -2ti:R 


27ir 
du  . dv 


_u 

r 


(3.66) 


-3T'  Yiz  = ^ + 


dz  di 


The  strain-displacement  matrix  [B]  for  triangular  element  and 
quadrilateral  element  are 


[B]  = [ B^  B2  Bj  ] 

[B]  = [ B,  B2  B3  B,,  ] 


where 
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Bi 


dN^ 

0 

1 

r 

dNi 

dz 


0 


dNj^ 

dz 

0 


dNj^ 

~dF 


1 !/■■■/  T2 


(3.67) 


From  equation  (3.63)  infinitesimal  area  in  plane  stress  is  dA 
= |J|dsdt  and  the  infinitesimal  volume  dv  may  be  similarly 
expressed  as  dv  = 27rr|J|dsdt.  The  element  stiffness  matrix 
[k]  of  the  axisymmetrically  loading  ring  element  can  be 
expressed  as 


[k]  = /^[B]‘"[D]  [B]dv  = [B]  ^ [D]  [B]  |j|rdsdt  ^2.68) 

Similarly,  the  consistent-mass  matrix  [m]  is 

[m]  = J p [N]  [N]  dv  = [N]  [N]  Ijjdsdt  (3.69) 

Also,  equivalent  nodal  loads  due  to  the  body  forces  are 


[fv]  = J^[N]'^[B^]dv  = 27ij'^^J^^[N]'^[B^]rdsdt  (3.70) 

Solution  of  the  Static  System 
The  solution  of  the  static  system 

[K]  {Ujj}  = F (3.71) 

is  an  important  step  in  the  finite  element  method  of  solution. 
The  Gaussian  elimination  method  is  used  to  solve  the  system 
equations.  Gaussian  elimination  can  be  reformulated  in  two 
phases,  factorization  and  substitution.  [L] [U]  decomposition 
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process  can  be  described  in  the  following  example.  A 
decomposition  algorithm  can  be  developed  to  obtain  all  the 
terms  of  [L]  and  [U]  from  the  original  terms  of  [K] . All 
terms  of  [L]  and  [U]  are  stored  in  the  original  matrix  [K]  as 
follows 


Uii 

a 

to 

• 

• 

L21 

U22 

• 

• 

• ^2n 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

.^1 

• 

• 

• 

^,n-i 

Un,n. 

(3.62) 


The  system  stiffness  matrix  can  be  decomposed  into  [L] [U]  as 


Ki2 

Ki3  • • 

■ 1 

0 

0 

• •* 

Uii 

U12 

Ul3  • • 

K21 

K22 

K23  • • 

L21 

1 

0 

• • 

0 

U22 

U23  • • 

K31 

K32 

K33  • • 

L31 

^32 

1 

• • 

0 

0 

U33  • • 

It  is  easy  to  construct  a row  of  [L]  and  a column  of  [U]  in 
succession  as  shown  in  table  3-1.  Such  an  algorithm  is  well 
adapted  to  a skyline  matrix.  The  operations  are  processed  in 
rows  for  [L]  and  in  columns  for  [U]  and  this  is  the  same 
sequence  when  terms  of  [L]  and  [U]  stored  in  the  skyline 
method  which  will  be  discussed  later. 
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Table  3-1  Construction  of  a row  of  [L]  and  a coluion  of  [U] 
in  succession 


index 

row 

in 

[L] 

column 
in  [U] 

from  which 

Lij  or  Uij  we 
get 

S=1 

1 

1 

Un=K” 

Uii=Kii 

s=2 

2 

1 

L2i=Kpi/Un 

1 

2 

Ui2=K,2 

2 

2 

LpiU,p+Upp=Kpp 

Upp=Kpp-Lp,Uip 

s=3 

3 

1 

In.=K,./S„ 

3 

2 

i^1^12'*‘i^2^22“K32 

^2“(K32~i^lUl2) 

/U22 

1 

3 

2 

3 

3 

3 

U 1 3+ L32U23+U33 

=K,3 

^33“K33-L3iUi3- 

^^2^23 

For  any  value  of  s 


Lsi  = 


i-l 

Ksi  - E ^smU„ 


/U„  i=l,2, . . . ,s-l 


(3.74) 


Ujs  = 

m=“l 

Matrices  [L]  and  [U]  are  stored  in 
and  the  previous  algorithm  becomes 


j=l,2 S (3.75) 

[K]  as  shown  in  Eq.(3.73) 


for  s=2 , 3 , . . . ,n 

for  i=l , 2 , . . . ,s-l 

for  m=l , 2 , . . . i-l 

K.;  = K . - K K . row  of  [L] 

SI  SI  sm  ml  L j 

Kic  = K.  - K-  K column  of  [U] 

IS  IS  im  Tns  >■  J 

normalize  rows  of  [L] 


Ksi  = K3.  / K,. 
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m=l , 2 


/ • • • 


s-1 


K, 


ss 


diagonal  term 


After  decomposition  by  the  above  algorithm,  the  matrix 
[K]  contains  all  the  terms  of  [L]  and  [U]  except  the  unit 
terms  on  the  diagonal  of  [L] . The  solution  of  the  system 
proceeds  in  two  steps  as  follows: 

(1)  Forward  substitution  in  the  lower  triangular  system. 
After  factorization  of  [K] , the  system  equation  may  be 
expressed  as  [L][U]{C)  = {F}.  The  solution  of  {U}  can  be 
obtained  first  by  solving  equation  [L]{G}={F).  The 
algorithm  is 

for  i=2 , 3 , . . .n 
c=0 

for  j=l , 2 , . . . , i-1 
C = C + K..F. 

F.  = F.  - C 

(2)  Backward  substitution  in  the  upper  triangular  system. 
After  {G}  is  solved,  the  system  displacement  may  be 
obtained  by  solving  [U]{C}={G}.  The  algorithm  is 


/ • • • / 


1 


c=0 


for  j=i+l 


, . . . , n 


c = c + KjjFj 
Fi  = (Frc)  / K-. 
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The  most  efficient  strategy  for  the  storage  of  global 
matrices  is  the  skyline  method.  It  consists  of  storing  the 
terms  of  [K]  or  [M]  by  rows  and  columns  of  variable  length. 
The  skyline  is  the  envelope  of  the  column  of  variable  heights. 
The  envelope  is  symmetrical  even  when  [K]  or  [M]  is  non- 
symmetrical.  The  sequence  of  [U]  in  the  columns  from  top  to 
diagonal  is  the  same  sequence  of  [L]  in  the  rows  from  left  to 
diagonal  for  same  s=i.  To  avoid  operations  by  zero  outside 
the  skyline  profile,  the  factorization  algorithm  must  be 
modified  as  follows 
for  s=2 , 3 , . . . ,n 

for  i=ios, . . . ,s-l 

for  m=max(i^.,i^J  , . . .i-1 


Ksi 

= Ksi 

- 

row  of 

[L] 

Kis 

= K,3 

- KiAs 

column 

Of  [U] 

Kg.  = Kgj  / K..  normalize  rows  of  [L] 

Kss  = Kgg  - Kg^K^g  diagonal  term 

where  i^.  and  i^g  are  the  integer  identification  of  the  high 
terms  in  columns  s and  i.  The  same  value  identifies  the  left- 
most terms  of  rows  i and  s.  The  terms  outside  the  i or  i . 

os  01 

are  zero  and  only  the  terms  within  i„^  and  s or  i . and  i 

should  be  processed.  When  processing  terms  of  [L]  and  [U] , we 

should  select  the  maximum  value  of  i^g  or  i^.  that  is  the 

shorter  in  row  width  or  column  height.  If  the  [K]  or  [M] 

matrix  is  stored  in  the  full  matrix,  the  choice  of  i or  i . 

os  01 
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is  not  important  and  it  only  requires  more  operations  on  the 
K-.  or  My  terms  with  zero.  But  when  we  compact  and  store  the 
matrix  by  the  skyline  method,  it  is  very  important  to  avoid 
any  unnecessary  zero  operations , because  any  unnecessary 
operations  will  process  the  terms  of  [K]  or  [M]  that  are  not 
stored  and  this  will  introduce  error.  The  solution  procedure 
should  also  be  modified  to  avoid  operations  outside  the 
3]^y3_in©  profile.  The  forward  substitution  equation  should  be 
modified  as  follows 
for  i=2 , 3 , . . .n 
c=0 

for  j~i(3s' ' * * 

c = c + KyF. 

F,  = F.  - C 

In  the  backward  substitution  the  algorithm  should  also  be 
modified  as 

= VK,.. 

for  1 — n“l ,n*2 , . . . , 1 
c=0 

for  j=i+l , . . . ,n 
if  i >=  3os 

C = C + KyFj 

F.  = (F.-c)  / Ky 

The  i is  the  identification  of  the  first  non-zero  term  in  j- 

•*  OS 

th  column.  The  procedure  involves  multiplication  between  the 
i-th  row  of  [U]  and  the  vector  {F).  There  are  variable  column 
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heights  between  (i+l)-th  column  and  n-th  column  and  we  have  to 
check  that  the  index  i is  below  the  first  non-zero  term  of  j- 
th  column. 

In  the  case  of  symmetric  systems,  the  algorithm  should  be 
modified  again  to  avoid  computation  for  the  lower  triangular 
terms  of  [K]  or  [M]  and  their  storage.  The  algorithm  is 
modified  as 

for  s=2 , 3 , . . . ,n 

for  i^s+l'  • • • 

for  m=max  ( ioi  / ios ) ' • • • ' 

Kie  = K.  - K.  K 
IS  IS  im  ^s 

c=0 

for  m=i^g, . . . ,s-l 

c c 

K.S  = 

Kss  = K33  - C 

Solution  of  the  Dynamic  System 
The  equation  of  equilibrium  governing  the  linear  dynamic 
response  of  a system  is  described  as 


[K]  {U}  + [C]  {U}  + [M]  {U}  = {R}  (3.76) 

The  procedure  used  in  this  study  to  solve  the  dynamic  system 
equation  is  the  direct  integration  method  [55,56].  In  the 
direct  integration  method  the  equation  (3.76)  is  integrated 
using  a numerical  step-by-step  procedure.  Instead  of 
satisfying  equation  (3.76)  at  any  time  t,  it  is  aimed  to 
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satisfy  equation  (3.76)  at  discrete  time  intervals  At  apart. 
This  means  that  equilibrium  is  sought  at  discrete  time  points 
within  the  intervals  of  solution.  In  some  cases  the  accuracy, 
stability,  and  cost  of  the  solution  procedure  is  critical  to 
the  time  interval  At.  It  is  very  important  to  select  one 
algorithm  with  good  characteristics.  The  Houbolt,  Wilson,  and 
Newmark  methods  use  the  equilibrium  conditions  at  time  t+At 
and  are  called  implicit  integration  methods.  They  are 
unconditionally  stable  integration  schemes.  The  effectness  of 
the  unconditionally  stable  integration  scheme  is  derived  from 
the  fact  that  in  order  to  obtain  accuracy  in  the  integration, 
the  time  step  At  can  be  selected  without  a requirement  for  a 
time  step  formula  of  critical  condition.  The  Houbolt 's  method 
is  unconditionally  stable,  but  it  is  not  self -starting.  Other 
methods  must  be  used  to  start  a solution.  The  Wilson  and 
Newmark  methods  are  also  unconditionally  stable.  In  this  study 
the  Wilson  method  is  selected  to  solve  the  dynamic  system 
equations.  The  Wilson  0 method  is  an  extension  of  the  linear 
acceleration  method.  In  the  Wilson  6 method  the  acceleration 
is  assumed  to  be  linear  from  time  t to  time  t+0t,  where  0 > 
1.0.  For  unconditional  stability  we  need  to  use  0 > 1.37  and 
the  optimum  value  0=1.420815  is  employed  in  the  program.  Let 
T denote  the  increase  in  time,  where  O<r<0At,  then  for  the 
time  interval  t to  t+0At,  it  is  assume  that 
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= {U},  -H 


(3.77) 


Integrating  equation  (3.77),  we  obtain 


{U},„  = {U)t  . {U)t  + ^^((U)c«4,-{U)t) 


(3.78) 


Integrating  equation  (3.78),  we  obtain 


Using  equations  (3.78)  and  (3.79),  we  have  at  time  T=0t 


(3.80) 


{U)t..At  = {U},+0At{U},+  ®^({U},,,^,+2{U},) 
from  which  we  can  solve 

To  obtain  the  solution  for  the  displacements,  velocities,  and 
accelerations  at  time  t+At,  the  equilibrium  equation  (3.76)  is 
considered  at  time  t+0At.  The  acceleration  is  assumed  to  vary 
linearly  and  a linearly  projected  load  vector  is  used.  The 
equation  (3.76)  becomes 
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• • 


[M]  {U}  t.0At 


+ [C]  {U} 


t+0At 


+ [K]  {U} 


where 


{R},,eAt  = 0({R}t.6At-tR>t) 


(3.85) 


Substituting  equations  (3.82)  and  (3.83)  into  equation  (3.84), 
we  can  obtain  an  equation  to  solve  the  displacement  at  time 
t+0At.  Then  substituting  the  result  into  equation  (3.68)  we 
obtain  the  acceleration  at  t+6At  which  is  used  in  equations 
(3.77),  (3.78),  and  (3.79)  to  calculate  accelerations, 
velocities,  and  displacements  at  time  t+At.  Then  the  process 
can  go  forward  for  the  next  time  step. 

The  next  important  factor  in  the  solution  of  the  dynamic 
system  is  choosing  the  time  interval  At.  From  an  efficiency 
study  of  implicit  and  explicit  integration  algorithms  by 
Katona,  Thompson,  and  Smith  [57],  a stability  criteria  can  be 
written  as 


where  c is  the  speed  of  sound  in  the  medium.  For  an  elastic 
solid  c is 


c 


E(l-u) 


(3.87) 


\ p (1+u)  (l-2u) 


with  E equal  to  Young's  modulus,  u equal  to  Poisson's  ratio, 
and  1 is  the  smallest  dimension  of  an  element.  For  the  Wilson 
method,  b is  equal  to  1/3. 
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Automatic  Processing  of  the  Boundary  Condition 

During  the  design  iteration,  the  shape  boundary  is 
updated  and  smoothed.  The  mesh  generation  has  to  be  executed 
for  each  iteration.  In  this  condition,  we  have  to  update  the 
constraints  and  loads  after  mesh  generation.  It  is  not  easy 
for  the  user  to  prepare  the  data  for  the  constraints  and  loads 
because  the  user  must  go  through  coordinates  and  connectivity 
tables  to  identify  elements  which  carry  loads  or  have 
constraints  on  it.  In  this  study  special  algorithms  are 
developed  to  handle  the  boundary  conditions  for  the  two- 
dimensional  finite  element  case  automatically.  Geometric 
elements  which  are  used  to  described  the  shape  boundary  are 
also  used  for  automatic  handling  of  the  boundary  conditions. 
The  physical  data  is  combined  with  geometric  data  to  represent 
a boundary  condition.  For  example,  the  point  load  can  be 
represented  by  combining  a geometric  data,  point  with 
coordinates  (x,y),  and  a physical  data,  the  loads  in  x and  y 
direction. 

In  the  design  process  the  boundary  condition  of  the 
finite  element  method  must  be  kept  the  same.  In  this  study  we 
employ  the  elimination  algorithm  and  this  changes  the  shape  of 
the  structure.  The  algorithms  for  the  automatic  handling  of 
the  boundary  condition  guarantee  the  same  boundary  conditions 
during  the  design  iteration.  Any  element  has  a node  or  a side 
located  on  the  constraint  element,  load  element,  or  prescribed 
displacement  element  will  not  be  removed  from  the  structure 
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When  the  elimination  algorithm  requires  it.  The  method 
developed  in  this  study  will  keep  the  physical  model  of  the 
two-dimensional  finite  element  method  consistent  during  the 
design  iterations  and  it  provides  a natural  and  easy  way  to 
represent  the  physical  model  of  finite  element  in  computer. 
Description  of  Constraints 

The  constraints  contain  the  nodal  degree  of  freedom  that 
is  fixed  in  the  finite  element  method.  In  this  study,  we  add 
one  requirement  for  the  constraints.  There  are  sometimes 
parts  of  the  structure  which  should  be  maintained  even  when 
the  load  carried  in  that  area  is  low.  The  elimination 

algorithm  will  remove  material  with  low  stress,  but  this  area 
is  not  fixed  and  will  be  processed  in  the  finite  element 
calculations.  Therefore,  we  define  a flag  with  the  geometric 
data  structure  of  the  constraint  element  to  handle  this 
problem.  The  data  structure  defining  constraint  elements  is 
listed  in  table  3-2. 


Table  3-2  Data  structure  of  the  constraint  representations 


type 

flag 

data 

data 

data 

data 

data 

flag 

flag 

point 

0 

X 

y 

0 

0 

0 

0, 

1 

0, 

1 

line 

0 

xl 

yi 

x2 

yi 

0 

0, 

1 

0, 

1 

arc 

1,-1 

xc 

yc 

r 

01 

02 

0, 

1 

0, 

1 

box 

0 

xl 

yi 

XU 

yu 

0 

0, 

1 

0, 

1 

The  first  seven  data  fields  are  defined  the  same  as  the 
geometry  element  which  describes  the  boundary  of  shape  and 
this  maintains  consistency  in  the  program.  The  physical  data 
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here  are  the  flags  used  to  indicate  constrained  degrees  of 
freedom  in  the  x-axis  or  y-axis.  If  we  only  want  to  keep 
material  in  structure,  the  flags  are  both  set  to  zero  for  the 
X and  the  y axis. 

After  the  mesh  generation  and  node  renumbering,  the 
program  must  prepare  the  data  for  the  finite  element  analysis 
and  the  boundary  conditions  should  be  updated  first.  But  the 
problem  here  is  to  obtain  constraint  element  identifications 
for  every  node.  The  algorithm  used  in  this  study  is  similar 
to  the  algorithm  used  in  interactive  two-dimensional  computer 
graphics  to  pick  a character,  a geometric  element,  or  a symbol 
on  the  display  screen.  The  coordinates  of  the  mesh  node  are 
used  as  the  world  coordinates  of  a locator  in  the  computer 
graphics.  We  must  decide  which  node  locates  on  which 
constraint  elements. 

The  bovinding  box  algorithm  is  used  for  the  same  reasons. 
When  we  want  to  decide  whether  a point  is  located  on  a 
geometry  element  or  not,  the  distance  between  node  and  the 
geometry  element  must  be  calculated.  But  it  is  costly  to 
calculate  the  distance  between  the  node  and  every  element  when 
the  number  of  elements  is  large.  We  can  use  a line  as  an 
example.  We  are  dealing  with  the  line  segment  instead  of  the 
line.  When  the  distance  between  the  node  and  the  line  is 
zero,  this  does  not  mean  that  the  node  is  located  on  the  line 
segment  as  shown  on  figure  3-3.  Further  checking  is  needed  to 
confirm  that  result.  If  we  go  through  the  distance 
calculation  and  further  check  for  every  line  segment,  the 
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computer  time  used  will  increase  rapidly  with  the  number  of 
line  segments  increased.  The  bounding  box  algorithm  will 
reduce  unnecessary  distance  calculation  for  impossible  line 
segments  and  simplify  the  conditions  to  check  location. 

The  bounding  box  of  a geometric  element  is  a rectangular 
box  that  can  contain  whole  element  and  the  distance  tolerance 
is  added  to  it  as  shown  in  figure  3-4.  If  the  bounding  box  of 
element  does  not  contain  the  node,  the  node  can  never  be 
located  on  this  element  and  it  is  not  necessary  to  calculate 
the  distance.  The  condition  for  the  bounding  box  containing 
a node  is  that 


x„.  < X < x„^„ 

mm  max 


V • < y < y 

•^min  — 1 — -I  max 


If  the  bounding  box  contains  the  node,  the  distance 
calculation  and  further  checks  then  take  place.  If  the  node 
locates  on  more  than  one  constraint  element,  the  constraint 
flags  of  elements  are  added.  If  the  value  of  the  flag  is 
greater  than  zero,  the  node  is  fixed  and  we  must  update  the 
system  matrix  of  finite  element.  Otherwise  the  degree  of 
freedom  of  node  is  not  constrained.  If  the  node  does  locate 
on  the  constraint  element,  but  the  constraint  flag  is  zero. 
In  this  condition,  the  finite  element  that  contains  the  node 
can't  be  eliminated  from  the  structure  and  there  is  no  need  to 
update  the  system  matrix. 

Representation  of  Loads 


The  point  load  is  defined  by  the  coordinates  of  a point 
(x,y)  and  the  loads  applied  on  the  node  in  the  x and  y axis 
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direction  as  shown  in  figure  3-5.  The  bounding  box  algorithm 
also  used  here  to  speed  up  the  process.  The  coordinates  of  a 
point  within  the  element  can  be  interpolated  by 


M = . 

VI 

The  point  load  can  also  be  interpolated  as 


(3.88) 


{P} 


NiPXi 

NiPVi 


(3.89) 


For  a point  load  on  an  element,  we  only  know  the  coordinates 
in  the  (x,y)  coordinate  system. 

In  order  to  evaluate  the  contribution  of  the  point  load 
to  the  nodal  points  of  the  element,  we  first  need  to  calculate 
the  corresponding  natural  coordinates  of  the  point  load  within 
element.  For  a triangular  element 

X = (1-s-t)  xl  + s x2  + t x3  = xl  + (x2-l)  s + (x3-xl)  t 
=>  x-xl  = (x2-xl)s  + (x3-xl)t 
Similarly,  we  can  obtain 

y-yl  = (y2-yl)s  + (y3-yl)t 
In  matrix  form,  the  above  equations  becomes 


x2-xl  x3-xl 

y2-yl  y3-yl. 

w 

(3.90) 


Then  we  can  evaluate  (s,t)  as 


1 y3-iy 
|det|  - (y2-yl) 


-(x3-xl)  ^-xll 
x2-xl  V“yiJ 


(3.91) 
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where  |det|  = (x2-xl) (y3-yl)  -(y2-yl ) (x3-xl ) 

For  the  applied  point  load,  we  know  the  (x,y)  coordinates  and 
the  natural  coordinates  (s,t)  of  the  point  load  on  the  element 
can  be  found  by  Eq.(3.91).  The  contribution  of  point  load  to 
all  nodes  of  the  triangular  element  can  then  be  calculated  by 
Px.  = N.(s,t)  Px,  Py.  = N.(s,t)  Py  i=l,2,3  (3.92) 

Similarly,  for  quadrilateral  element,  we  can  use  equation 
(3.92)  to  get 

4x  = Ax  + Bx  s + Cx  t + Dx  s t 
4y  = Ay  + By  s + Cy  t + Dy  s t 
where  Ax=xl+x2+x3+x4 , Bx=x2+x3-xl-x4 , Cx=x3+x4-xl-x2 , 
Dx=xl+x3-x2-x4 

Ay=yl+y2+y3+y4 , By=y2+y3-yl-y4 , Cy=y3+y4-yl-y2 , 
Dy=yl+y3-y2-y4 

Using  the  niimerical  method,  we  can  solve  for  the  t value  first 
and  s value  can  be  obtained  by  using  the  above  equations.  In 
this  way,  we  can  automatically  add  point  loads  to  the  system 
equation  in  finite  element  formulation. 

The  surface  loading  element  in  this  study  are  assumed  to 
be  linear  and  to  act  on  the  only  the  bovindary  of  the  element. 
The  load  at  any  point  on  the  surface  loading  element  can  be 
calculated  by  linear  interpolation.  The  classifications  of 
pressure  elements  are  listed  in  the  table  3-3. 

The  geometric  data  use  he  same  classification  for  the 
geometric  elements  of  the  shape  boundary.  The  surface  loads 
at  the  two  end  points  are  used  to  interpolate  the  load  on  the 
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Table  3-3  Data  structure  of  surface  loads 


type 

flag 

dl 

d2 

d3 

d4 

d5 

Tx 

Ty 

Tx 

Ty 

line 

0 

xl 

yi 

x2 

y2 

0 

txl 

tyl 

tx2 

ty2 

arc 

-1,1 

xc 

yc 

r 

ei 

02 

txl 

tyl 

tx2 

ty2 

surface  loading  element.  The  bounding  box  algorithm  is  again 
used  to  check  if  the  side  of  the  finite  element  locates  on  any 
surface  loading  element  or  not.  The  equivalent  nodal  loads 
for  the  element  are  derived  in  the  following. 

Assuming  the  surface  load  acts  on  the  side  of  the  element 
as  shown  in  figure  3-6,  the  load  will  vary  linearly  from  point 
1 to  point  2.  Shape  function  on  the  side  has  a constant  value 
of  s or  t and  only  the  node  on  side  carrys  loads.  The  node 
not  on  this  side  has  no  load  at  all,  we  can  use  the  technique 
for  one-dimensional  problems  to  solve  the  equivalent  nodal 
forces  acting  on  the  two  end  points  of  the  side.  We  can 
therefore  assume  the  surface  loads  function  to  be 

T(s)  = CO  + Cl  s (3.93) 
The  boundary  condition  at  two  the  end  points  are 


T(0) 

= CO 

+ Cl 

0 

= T1 

T(l) 

= CO 

+ Cl 

1 

= T2 

We  can  then  obtain  the  coefficients  as  C1=T1,  C2=(T2-T1 )/l . 
The  linear  displacement  function  are: 

Nl(s)  = N2(s)=y  (3.94) 


The  equivalent  nodal  force  at  nodes  are 
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F,  = r^Nl  (s)T(s)ds  = (T1+  '^^■'^^s)ds  (3.95) 

Jo  Jo  L 1 

F,  = f^N2  (s)T(s)ds  = (T1+  '^^~'^^s)ds  (3.96) 

^ Jo  Jo  L 1 

We  can  evaluate  F^  and  F2  as 

F,  = — [3L(T1+T2) -1  (T1+2T2)  ] (3.97) 

6L 


F,  = 


6L 


[ ITI  + 21T2  ] 


(3.98) 


■'V 

II 

1 

1 

to 

6L 

3L-1  3L-21 
1 21 


a 


We  can  express  the  above  equation  in  matrix  form  as 

(3.99) 

For  the  case  of  L=1  and  constant  element  thickness,  we  have 

(3.100) 


[■ 

F, 


= h^ 


’2  il  fril 
1 2J  V2J 


The  procedure  checks  the  equivalent  nodal  force  and  updates 
the  system  force  vector  element  by  element.  The  program 
developed  here  can  handle  boundary  conditions  and  loads 
automatically  so  that  the  optimum  design  iterations  can  go 
from  beginning  to  the  end  without  user  interference. 

Postprocessor  for  the  Finite  Element  Method  Analysis 
In  order  to  quickly  evaluate  the  results,  a postprocessor 
is  also  developed  to  display  the  data  for  the  statics  or 
dynamics  analysis  on  the  screen.  The  postprocessor  can 
display  displacements,  stress  contour,  and 


stress 


stress 
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distribution.  The  display  of  velocity  and  acceleration  fields 
is  an  option  in  the  dynamic  analysis  program.  The  algorithm 
for  stress  contour  and  stress  distribution  uses  interpolation 
to  locate  the  contour  lines  which  pass  through  the  finite 
element.  First  we  have  to  define  the  stress  gradient  for  the 
contour.  If  there  is  at  least  one  stress  contour  that  pass 
through  the  element,  the  points  on  the  element  sides  should  be 
located  by  interpolating  the  stress  values  of  the  end  points 
of  the  side.  Then  the  element  is  split  by  connecting  those 
intersection  points  on  the  sides  and  the  process  is  repeated 
for  the  remaining  polygon. 
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Figure  3-1  Reference  element  and  real  element  for 


three-node  isoparameteric  triangular 
element 


Figure  3-2  Reference  element  and  real  element  for 


four-node  isoparameteric  quadrilateral 
element 
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Figure  3-3  Check  point  location  with  the  boundary 
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Figure  3-5  Point  load  acting  on  triangular  element 


Figure  3-6  Linear  surface  load  acting  on  the 
boundary 
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stress  contour  of  postprocessor 


CHAPTER  4 

SHAPE  OPTIMIZATION  AND  RULES 
Introduction 

The  elimination  concept  is  used  in  this  study  to  optimize 
the  structure  with  minimum  weight  under  static  and  dynamic 
loads  and  to  approach  a viniform  stress  distribution  within  the 
structure  as  closely  as  possible.  In  the  mathematical 
programming  approach,  the  control  points  on  the  parameterized 
boundaries  are  chosen  as  the  decision  variables  to  change  the 
shape.  For  complex  problems  where  the  shape  is  lanpredictable , 
it  is  not  easy  to  define  appropriately  parameters  for  the 
boundaries  of  shape.  Furthermore,  using  parameterized 
boundaries  to  optimize  shape  doesn't  allow  the  creation  of 
holes  inside  the  design  domain.  This  may  eliminate  a better 
alternative  design,  if  it  exists  by  constraining  the  movement 
of  control  points.  It  is  possible  that  the  best  optimum 
design  has  holes  inside  the  design  domain. 

A structural  shape  can  be  represented  completely  by  the 
elements  which  define  the  entire  domain  of  shape.  The  change 
of  shape  can  be  accomplished  by  eliminating  or  adding 
elements.  In  this  study  the  finite  element  method  is  coupled 
with  the  optimization  criterions  and  the  design  rules  to 
change  the  shape.  The  elements  which  represent  the  design 
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domain  are  also  the  elements  used  in  the  finite  element 
analysis.  The  results  of  the  finite  element  analysis  are 
naturally  utilized  to  change  the  shape.  Every  element  without 
constraints  can  be  removed  from  the  structure  and  there  is  no 
restriction  about  removing  elements  from  the  boundary  or  from 
the  interior  of  the  design  domain.  The  optimum  design  with 
holes  inside  structure  which  can't  be  achieved  by  methods 
which  parameterize  the  boundary  of  shape  will  be  naturally 
created  by  using  the  elimination  algorithm.  In  the  process  of 
elimination,  the  meshes  to  be  removed  are  grouped  into 
subdomains  of  shape  and  the  Boolean  set  operation  is  used  to 
update  the  shape. 

Elimination  Criteria 

For  the  case  of  static  or  inertial  loads,  the  results  of 
the  finite  element  analysis  are  used  directly  in  the 
optimization  algorithm.  In  the  impulsive  loading  cases,  the 
stress  at  the  of  same  point  is  different  at  different  time 
steps.  The  material  fails  when  the  dynamic  stress  exceeds  the 
strength  and  is  not  dependent  on  the  time  steps.  For 
impulsive  loads,  the  maximum  stress  in  an  element  at  any  time 
step  is  saved  for  the  duration  response  history.  In  this 
study,  it  is  assumed  that  no  damping  exists  in  the  system.  One 
cycle  is  defined  as  the  time  needed  for  a stress  wave  to 
travel  twice  the  distance  between  two  furthest  points  of  the 
structure.  The  response  period  used  to  determine  the  maximum 
nodal  stress  is  selected  to  be  four  cycles.  The  maximum  nodal 
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stresses  during  this  period  are  then  used  in  the  optimization 
process . 

The  characteristics  of  the  stress  distribution  at  every 
design  iteration  can  be  obtained  by  a histogram.  The  element 
used  in  the  finite  element  analysis  is  the  basic  unit  for  the 
elimination  algorithm.  After  the  finite  element  calculation, 
we  can  determine  the  stress  at  every  of  node.  The  same  node 
shared  by  the  adjacent  elements  may  have  different  stress 
values  and  the  arithmetic  average  is  taken  to  represent  the 
stress  at  the  node.  The  distortion  energy  theory  is  used  to 
predict  state  of  combined  stress  in  the  structure. 

Using  the  nodal  stresses  information  the  Mises-Hencky 
equivalent  stress  is  calculated  and  the  element  stresses  are 
grouped  into  different  stress  levels.  This  grouped  data 
constructs  the  frequency  distribution  of  the  element  stresses 
as  shown  in  figure  4-1.  Some  characteristics  of  location  and 
dispersion  can  be  used  to  appraise  the  stress  distribution 
within  design  domain. 

The  elimination  algorithm  removes  material  from  the 
structure,  but  there  are  some  areas  of  material  that  must  not 
be  removed  in  order  to  keep  the  physical  model  and  boundary 
conditions  consistent  with  original  status.  Any  element  that 
can't  be  removed  from  the  structure  is  called  a constrained 
element  and  they  fall  into  one  of  the  following 
classifications : 
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(1)  The  element  carrying  loads.  The  element  that  has 
point  load  or  surface  pressure  on  it  must  be  kept  in  the 
structure . 

(2)  The  element  which  has  constraints  imposed  on  it.  The 
boundary  conditions  of  the  finite  element  model  must  also 
be  kept  without  change  during  the  design  iteration  in 
order  not  to  violate  the  boundary  constraints. 

(3)  User  designated  element.  The  elimination  algorithm 
progresses  naturally  to  remove  the  elements  carrying  low 
stress.  Sometimes  this  will  create  holes  inside  the 
design  domain  and  may  cause  manufacturing  difficulties. 
If  the  user  prefers  to  keep  some  part  of  the  material 
within  the  structure  without  removing  it,  a designated 
subdomain  of  structure,  a line,  a arc  or  a box,  can  be 
employed  in  the  program  to  keep  material  in  that 
subdomain.  Any  element  that  has  intersection  with  those 
designated  elements  can  not  be  eliminated. 

Any  one  of  the  above  three  elements  must  be  kept  in  the 
structure  no  matter  how  low  the  stress  it  carries.  When  we 
group  elements  into  different  stress  levels,  the  constrained 
elements  are  excluded  from  the  grouping  in  order  not  to 
distort  the  criterion  for  the  elimination.  If  we  count  in  the 
stress  values  of  the  constrained  elements,  they  will  influence 
the  frequency  distribution  and  this  effect  will  increase  with 
the  number  of  constrained  elements.  In  the  following 
discussion  about  the  grouping  and  elimination,  the  element 
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without  any  coininent  is  the  free  element  that  doesn't  belong  to 
any  of  the  above  three  categories. 

The  range  of  the  distribution  is  the  difference  between 
the  maximum  and  minimum  element  stresses.  In  order  to 
classify  the  data  into  some  categories  of  frequency 
distribution,  the  range  is  divided  into  intervals.  The 
intervals  represent  the  different  stress  levels  for  placing 
elements  into  categories.  The  number  of  intervals  is  equal  to 
the  number  of  categories  and  the  selection  of  this  number 
affects  the  process  of  elimination.  If  a small  category 
number  is  used,  the  number  of  elements  in  each  category  will 
be  quite  large.  The  speed  of  elimination  may  be  fast,  but  the 
discontinuity  condition  may  occur  frequently  during  the 
process.  Although  the  program  can  handle  the  discontinuity 
condition,  this  is  not  a desirable  condition  to  happen  during 
the  optimization  process.  On  the  other  hand,  if  the  number  of 
categories  is  large,  there  will  only  be  few  elements  in  each 
category  and  the  optimization  process  will  become  very  slow. 
The  category  number  used  in  the  program  is  thirty.  The  range 
is  then  divided  into  thirty  intervals.  The  elements  are 
accordingly  classified  into  these  categories  and  the  frequency 
distribution  of  element  stresses  is  therefore  established  for 
further  analysis. 

When  we  extract  data  from  the  developed  histogram  the 
statistical  formula  [58]  for  the  grouped  data  should  be 
applied  where  every  observation  is  set  equal  to  its  category 
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mark.  The  category  mark  used  here  is  the  average  stress  of 
elements  within  the  same  category.  The  sample  mean  of  N 
categories  is 


where  x.  is  the  i-th  category  mark  and  f ■ is  the  corresponding 
frequency.  Measures  of  the  center  location  describe  one 
important  aspect  of  the  data,  their  average,  but  they  provide 
no  information  about  other  basic  characteristics.  We  require 
ways  of  measuring  the  dispersion  of  data  and  the  statistical 
measures  which  provide  this  information  are  called  measures  of 
variation.  The  dispersion  of  a set  of  data  is  small  if  the 
values  are  closely  bunched  about  their  mean  and  is  large  if 
the  values  are  scattered  widely  about  their  mean.  It  would 
seem  reasonable  to  measure  the  variation  of  a set  of  data  in 
terms  of  the  amounts  by  which  the  values  deviate  from  their 
mean.  The  variance  of  a sample  is  usually  defined  as  the  mean 
of  the  squared  deviations  from  the  mean  as 


N N 


(4.1) 


N 


(4.2) 


a 


2 - i=l 


n-1 


The  standard  deviation  of  N categories  with  mean  is  given  by 
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E (Xi-|i)2f.  (4.3) 
i=l 

n-1 

From  the  definition  of  the  standard  deviation  that  if  the 
standard  deviation  of  a set  of  data  is  small,  the  values  are 
concentrated  near  the  mean.  If  standard  deviation  is  large, 
the  values  are  scattered  widely  about  the  mean.  For  an  ideal 
condition  of  constant  stress  distribution,  the  range  and 
standard  deviation  are  both  very  small  in  the  histogram.  The 
is  the  goal  of  the  elimination  process  is,  therefore,  to 
approach  a good  frequency  distribution  in  the  histogram. 

The  stress  criterion  used  to  eliminate  elements  is  the 
category  mark  that  has  the  largest  increase  in  frequency  or 
number  of  elements  for  two  consecutive  categories.  This  also 
indicates  that  there  are  a number  of  elements  having  almost 
the  same  stress  level  in  the  following  category.  If  we 
eliminate  too  many  elements  from  the  structure,  the  first 
drawback  is  that  a discontinuity  may  likely  happen  in  the 
process  of  changing  shape  and  it  is  more  serious  when  the 
elements  belong  to  one  of  the  major  groups  that  carry  loads  in 
the  structure.  When  the  following  group  has  a significant 
number  of  elements  in  it,  then  the  group  plays  an  important 
role  in  transmitting  loads  to  a good  percentage  of  the 
structure.  If  we  remove  the  following  group,  the  stress 
distribution  will  increase  substantially  because  it  is  one  of 
the  important  groups  since  it  occupies  a good  percentage  of 
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the  volume  and  carries  higher  loads.  The  elements  always  tend 
to  concentrate  near  the  mean  and  it  is  not  unusual  to  obtain 
a large  increase  in  number  of  elements  around  the  mean.  In 
order  to  avoid  removing  major  stress  distribution  bands  in  the 
structure,  we  seek  the  largest  increase  in  element  number  of 
two  consecutive  categories  after  one  standard  deviation.  All 
elements  with  stress  below  this  category  will  be  eliminated 
from  the  structure. 

Elimination  Groups 

After  we  determine  the  element  identifications  for  the 
elimination,  those  element  may  be  randomly  located  inside  the 
design  domain.  In  order  to  remove  an  element  from  the  design 
domain,  we  need  a union  Boolean  operator  between  structure  and 
the  eliminated  elements.  If  we  remove  material  from  the 
structure  element  by  element  for  all  candidate  elements,  more 
computer  time  will  be  needed  and  the  conditions  for  Boolean 
operator  also  become  more  complex.  A special  Algorithm  is 
developed  for  the  elimination  of  elements  which  are  grouped  by 
having  common  sides  with  other  candidate  elements.  The  result 
of  this  grouping  technique  converts  a large  number  of  elements 
marked  for  elimination  into  a few  elimination  subdomains. 

We  can  define  a two-dimensional  objects  by  the  boundaries 
of  the  objects  which  are  constructed  by  geometric  elements 
linked  in  a counterclockwise  direction.  This  is  a boundary 
representation  model  illustrated  in  figure  4-2.  In  general 
case,  the  shape  of  the  structure  which  is  defined  as  the 
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design  domain  is  a combination  of  solid  areas  and  holes.  The 
elimination  group  is  defined  as  an  object  that  is  an  empty 
object  or  a hole.  In  the  optimization  iterations,  the  shape 
of  the  previous  iteration  and  the  objects  of  the  current 
elimination  groups  are  combined  by  a union  Boolean  set 
operator  to  create  the  optimum  shape  for  the  current 
iteration . 

A group  of  elements  may  be  regarded  as  a number  of 
elements  that  are  jointed  together  by  a common  side  between 
the  elements  and  they  are  enclosed  by  a closed  boundary  as 
shown  in  figure  4-3.  The  identifications  of  the  elimination 
elements  as  determined  by  the  optimization  criterion  are 
generally  arranged  in  a random  way.  There  is  no  information 
about  the  geometric  relationship  between  elements.  In  order 
to  group  elimination  elements,  we  start  with  s single 
elimination  element.  The  first  elimination  group  contains 
only  one  element  at  the  beginning.  Then  the  procedure  check 
the  remaining  elimination  elements  to  group  them  by  the 
following  rules: 

(1)  If  the  elimination  element  has  no  common  side  with 
any  existing  elimination  groups,  it  is  regarded  as  the 
first  one  in  a new  elimination  group. 

(2)  If  the  elimination  element  has  at  least  one  common 
side  with  an  existing  elimination  group,  it  is  jointed 
with  that  group  which  becomes  a new  and  updated  group. 

(3)  If  the  elimination  element  has  common  sides  with  more 
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than  one  of  the  existing  groups,  all  adjacent  groups  are 
jointed  together  by  this  elimination  element  and  this 
become  one  group.  The  counter  and  other  group  data 
structure  also  need  to  be  updated. 

After  this  procedure,  the  elimination  elements  are  divided 
into  groups  which  are  geometrically  jointed  together.  In 
order  to  define  them  as  hole  objects,  we  need  to  define  the 
boundary  which  encloses  the  elimination  elements. 

As  we  mentioned  earlier  the  object  is  defined  by  the 
boundary  in  a counterclockwise  seguence,  the  nodal  seguence  of 
the  finite  elements  is  also  defined  in  a counterclockwise 
direction  and  this  helps  us  to  trace  the  outline  of  the 
elimination  group.  A start  point  is  therefore  needed  to  begin 
the  tracing  procedure.  The  reguirement  of  this  start  is  that 
it  must  be  located  on  the  boundary  of  the  group.  The 
procedure  for  searching  for  a start  point  is  to  find  an 
extreme  point  which  is  located  on  boundary  and  is  used  by  one 
element  only.  We  search  for  the  point  with  minimum  x 
coordinate  first,  then  the  point  with  minimum  y coordinate  is 
sorted  out  to  find  the  start  point.  Any  point  on  the  boundary 
can  be  only  used  by  two  sides  of  adjacent  elements  when  both 
sides  are  also  located  on  the  boundary.  In  order  to  trace  the 
boundary,  we  must  find  out  all  elements  inside  the  elimination 
group  that  have  same  common  point  which  we  call  a tracing 
point,  like  the  start  point.  The  boundary  of  the  elimination 
group  is  accordingly  created  by  the  following  rules: 
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(1)  If  there  is  only  one  element  using  the  tracing  point, 
the  next  nodal  point  of  the  element,  in  the 
counterclockwise  sequence,  becomes  the  new  tracing  point. 

(2)  When  there  are  more  than  one  elements  using  the 
current  tracing  point,  all  the  next  nodal  points  of  those 
elements  are  candidates  for  the  tracing  point.  From 
those  candidates,  we  choose  the  tracing  point  by  checking 
the  requirement  that  the  current  and  next  tracing  point 
must  be  located  on  the  boundary.  This  means  that  the 
side  constructed  by  the  current  and  next  tracing  point  is 
used  by  one  element  only.  From  this  requirement  we  can 
find  out  the  next  tracing  point. 

The  procedure  is  repeated  until  the  start  point  becomes  the 
next  tracing  point  again,  as  shown  in  figure  4-3.  After  we 
obtain  all  the  boundaries  for  the  elimination  group.  It  will 
be  ready  for  the  union  Boolean  operator  to  create  the  new 
shape  by  which  the  elimination  groups  are  removed  from  the 
structure . 

Boolean  Operations  Between  Elimination  Groups  and  Shape 
At  the  beginning  of  the  design  iteration,  the  shape  of 
the  structure  is  a two-dimensional  object  with  or  without 
holes  inside  of  it.  In  the  boundary  representation  model,  the 
two-dimensional  objects  are  defined  by  the  boundaries  which 
are  constructed  by  the  geometric  elements  linked  in  a 
counterclockwise  sequence.  When  we  try  to  perform  the  Boolean 
operation  between  two  objects,  it  is  necessary  to  traverse  the 
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boundaries  of  the  objects  in  both  directions.  The  information 
of  the  boundaries  of  object  are  defined  by  doubly  linked  lists 
and  a special  flag  is  used  as  a solid  or  hole  indicator  of  the 
boundary  loops.  Since  the  boundary  of  a solid  or  a hole  is 
always  arranged  in  a counterclockwise  sequence,  we  need  a 
indicator  to  separate  a solid  and  a hole.  A Boolean  set 
operations  algorithm  that  can  be  used  to  unite,  intersect,  or 
subtract  objects  is  an  essential  part  of  updating  the  shape 
after  elimination.  Set  operations  offer  a tool  for  describing 
complex  objects  in  terms  of  a series  of  operations  on  simpler 
components.  In  this  study,  the  union  set  operation  is  applied 
to  change  the  shape  of  the  structure  by  eliminating  materials 
and  the  operations  are  used  mainly  in  uniting  the  boundaries 
of,  the  objects. 

The  data  structure  used  to  represent  a boundary  of  an 
object  is  a doubly  linked  list  which  describes  every  node  on 
the  boundary  in  a counterclockwise  sequence.  The  input  of  the 
Boolean  set  operation  algorithm  will  consist  of  two  doubly 
linked  lists,  A and  B.  Given  two  boundaries  of  objects  A and 
B,  the  collection  of  four  segments  AinB,  AoutB,  BinA,  and 
BoutA,  as  shown  in  figure  4-4,  formed  by  splitting  A against 
B and  symmetrically  B against  A are  called  the  boundary 
classification  of  A and  B.  From  the  boxindary  classification 
of  A and  B,  the  Boolean  combination  is  computed  as; 

A u B = AoutB  ® BoutA 

Assuming  A is  a solid  object  and  B is  a hole,  the  union 
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Boolean  set  operation  can  be  illustrated  by  the  figure  (4.)* 
In  the  figure,  the  segment  AinB  and  BinA  are  not  displayed 
because  B is  a hole  and  it  is  physically  empty  there. 

In  the  program,  the  intersection  segments,  AinB  and  BinA, 
are  calculated  first.  Segments  AoutB  and  BoutA  can  then  be 
obtained.  Then  we  apply  union  Boolean  set  operation  to  unite 
A and  B and  the  result  becomes  the  new  boundary  of  the  object. 
The  elimination  group  may  be  regarded  as  a hole  as  shown  in 
figure  4-5.  The  shape  of  the  structure  is  united  with  every 
elimination  group  to  form  the  new  shape,  as  shown  in  figure  4- 
5.  Because  the  shape  of  structure  is  a combination  of  solids 
and  holes,  the  object  defined  by  the  elimination  group  must  be 
checked  for  intersection  with  every  boundary  of  structure.  If 
the  elimination  object  has  more  than  one  intersection  with  the 
boundaries  of  structure,  the  union  operation  must  be  applied 
to  those  boundaries  one  by  one.  If  the  elimination  object 
does  not  intersect  with  any  boundary  of  the  shape,  it  is 
regarded  as  a new  hole  and  becomes  a part  of  the  structure. 

Discontinuity  Condition  in  Elimination  Optimization 

Material  discontinuity  within  the  structure  which  may 
occur  in  the  elimination  process  is  not  allowed  in  the  shape 
optimization.  The  elements  to  be  deleted  are  grouped  into 
different  stress  level . When  we  apply  the  Boolean  set 
operation  after  grouping  elements , discontinuity  may  occur 
during  the  process.  In  order  to  satisfy  the  continuity 
requirement,  a rule-based  algorithm  is  developed  to  handle 
this  problem. 
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The  basic  logic  used  in  detecting  the  discontinuity 
condition  is  that  the  boundaries  of  the  eliminated  groups  have 
at  least  two  intersection  segments  with  the  same  shape  loop. 
The  discontinuity  condition  with  three  or  more  intersection 
segments  is  not  dealt  with  in  this  study.  Assuming  that  a 
eliminated  group  has  two  intersection  segments  with  the  shape 
loop,  as  shown  in  figure  4-6,  the  discontinuity  will  occur  if 
we  remove  the  entire  group  from  the  structure. 

If  the  region  on  the  structure  is  small  and  there  is  no 
constraints  imposed  on  it,  we  may  remove  this  region 
simultaneously  with  the  eliminated  group  from  the  structure  as 
illustrated  by  the  dashed  line  region  in  figure  4-7.  In  such 
case  the  structure  is  divided  into  two  regions  by  the 
eliminated  group.  The  number  of  elements  in  those  two  regions 
may  not  be  easy  to  determine  and  we  need  only  to  check  the 
number  of  elements  on  the  boundary  for  each  region.  If  the 
region  has  only  few  elements  on  the  boundary  and  there  is  no 
constraints  on  this  region,  we  can  then  remove  all  the 
material  in  this  region  from  the  structure. 

If  there  are  some  constraints  on  both  divided  regions,  we 
allow  the  eliminated  group  to  divide  the  structure  into  two 
substructures.  We  must  keep  part  of  the  elements  in  the 
eliminated  group  to  connect  these  two  regions.  The  material 
to  be  kept  may  be  left  in  the  center  of  elimination  group  or 
near  the  intersection  boundaries.  In  this  study  we  use 
central  elements  of  the  eliminated  groups  to  allow  some 
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material  to  connect  the  two  regions  as  shown  in  figure  4-8. 
The  width  of  the  connecting  central  elements  should  be  at 
least  one  element  in  size.  This  size  is  not  always  the  same 
and  the  width  of  the  central  part  depends  on  the  width  of  the 
eliminated  group. 

Updating  the  Shape  of  the  Structure 
The  boundaries  of  the  shape  are  a combination  of  line, 
arc  and  Bezier  spline.  In  order  to  smooth  the  zigzag  boundary 
after  elimination,  the  Bezier  spline  is  applied  to  all  the 
zigzag  boimdaries  which  do  not  has  any  constraints  or  loads  on 
it.  The  shape  of  the  previous  iteration  is  used  as  the 
reference  for  updating  the  new  shape.  Any  point  that  is 
located  on  the  line  or  arc  in  the  previous  iteration  will  be 
classified  into  a line  or  a arc  according  to  the  type  of  the 
geometric  element  on  which  the  point  is  located.  All  those 
points  are  masked  by  a special  flag  which  points  to  the  data 
of  the  geometric  element.  Any  other  point  which  is  not 
located  on  constraints  or  loads  is  masked  as  a free  point. 
Using  those  masked  information,  we  can  travel  new  shape 
boundaries  in  a counterclockwise  to  extract  the  geometric 
elements  from  boundary  and  update  the  data  structure  of  the 
boundary  of  the  shape.  These  connected  boundaries  become  the 
new  shape  of  structure  for  next  design  iteration  as  shown  in 
figure  4-9. 
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Figure  4-1  Histogram  for  elimination 


Figure  4-2  Boundary  representation  of  structure 
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Figure  4-3  Bovmdary  representation  of  eliminated 
group 


Figure  4-4  Boolean  operation  of  two-dimensional 
shape 
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Figure  4-5  Boolean  operation  between  structure  anc 
eliminated  group 


Figure  4-6  Discontinuity  created  by  eliminated 
group 
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Figure  4-7  Unconstrained  structure  removed  with 
the  eliminated  group 
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Figure  4-9  Updating  the  boundary  of  structure 

using  lines,  arcs,  and  Bezier  splines 


CHAPTER  5 

OPTIMIZATION  PROGRAM  STRUCTURE 
Introduction 

We  will  discuss  in  this  chapter  the  entire  design 
procedure  for  shape  optimization  by  the  elimination  algorithm. 
In  order  to  provide  the  user  with  easy  interaction  and 
flexibility,  there  are  few  options  in  the  program  to  choose 
from  according  to  the  specific  needs.  The  program  can 
continue  from  the  first  iteration  to  the  end  without  any 
interruption.  We  can  also  control  the  procedure  to  be 

executed  iteration  by  iteration  in  order  to  evaluate  each 
result.  If  the  result  is  not  satisfactory,  we  may  interrupt 
the  program  execution  and  make  some  modifications.  In  the 
dynamic  analysis,  the  time  increment  is  relatively  small  and 
it  may  take  a long  time  for  just  one  design  iteration.  If  we 
have  to  repeat  the  design  process  from  the  first  iteration, 
considerable  time  will  be  wasted.  In  order  to  avoid  this,  the 
program  stores  the  necessary  data  after  every  design  iteration 
and  the  user  can  restart  the  optimization  process  again  from 
any  iteration  number.  This  can  save  considerable  time  and 
speeds  the  shape  optimization  process. 

This  study  investigates  the  feasibility  of  an 
algorithmatic  approach  to  the  two-dimensional  shape 
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optimization  of  structures  subjected  to  static  or  dynamic 
loads.  Combining  the  finite  element  method  and  rule-based 
algorithms  for  elimination,  boundary  smoothing,  and  shape 
update  provides  a seguential  search  for  the  optimal  shape  of 
a mechanical  component  or  a structural  element.  The  problem 
to  be  treated  in  this  study  can  be  stated  as  follows: 

Given:  An  initial  two-dimensional  system  with  simple 
(primitive)  geometry,  the  applied  loads,  and  the 
constraints. 

Find:  The  shape  of  the  system  with  minimum  volume  and 
the  best  stress  distribution. 

Subject  to:  The  appropriate  behavioral  constraints 
without  violating  the  maximum  allowable  stress  or 
violating  the  continuity  of  the  material. 

The  developed  algorithms  may  be  classified  as  a rule- 
based  optimization  technigue  rather  than  a formal  mathematical 
programming  optimization.  In  most  of  the  literature  on  shape 
optimization  algorithms,  the  configuration  of  the  structure 
for  the  optimization  process  is  always  defined  before  the 
execution  of  the  algorithm.  Some  of  the  geometric  parameters 
are  chosen  as  the  design  variables  in  order  to  change  the 
shape  of  structure.  In  the  optimization  process,  a good 
starting  point  is  always  an  important  factor  for  the  success 
of  the  process.  The  initial  configuration  used  in  most  of  the 
studies  utilizing  such  approach  is  selected  based  on  the 
author's  experiences.  In  other  words  the  shape  is 
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preselected,  only  its  parameters  are  optimized.  If  we  want  to 
achieve  the  ultimate  goal  of  an  automated  shape  design,  it  is 
desirable  to  allow  the  configuration  to  evolve  in  a systematic 
way  automatically.  In  the  rule-based  algorithm  developed  in 
this  study,  the  only  constraints  on  the  shape  are  the  physical 
requirements  and  there  is  no  constraints  on  what  the  final 
shape  looks  like  before  and  during  optimization  process.  The 
program  still  supports  the  user-defined  constraints  to 
maintain  certain  segments  in  the  structure  regardless  of  the 
stress  level  they  are  subjected  to.  This  user-defined 
constraint  option  can  generate  alternative  solutions  to  fit 
the  basic  requirements  of  the  user.  The  requirement  may  be  to 
simplify  the  manufacturing  of  the  part  or  it  may  be  just  to 
force  another  design  option  based  on  the  user's  experiences. 
The  user  can  then  select  the  best  design  from  the  obtained 
results  with  or  without  the  user-defined  constraints. 

Preparation  of  Data 
Initial  Configuration  of  Design 

The  initial  configuration  needed  to  start  the  program  is 
a simple  geometry  which  is  judged  to  be  large  enough  to 
encampus  the  desired  part.  This  simple  geometry  may  be  the 
working  space  for  the  mechanical  element  in  the  machine,  the 
space  for  setting  the  structure,  or  the  simplest  shape 
satisfying  the  requirements  of  continuity  and  boundary 
conditions.  It  is  not  necessary  to  describe  the  shape 
precisely  or  have  a preconceived  notion  on  what  it  should  look 
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like.  The  algorithm  will  automatically  generate  the  shape. 
If  it  is  necessary  to  create  the  boundary,  the  algorithm  will 
define  it  naturally.  This  is  quite  different  from  the 
mathematical  programming  approach  which  is  based  on 
parameterized  boundaries. 

In  the  optimization  approaches  using  parameterized 
boundaries,  the  initial  configuration  must  be  carefully 
selected,  because  the  process  can  only  move  the  control  points 
on  the  parameterized  boundaries  and  there  is  no  potential  for 
automatically  generating  the  new  parameterized  boundaries. 
The  result  will  be  closely  related  to  the  initial 
configuration  and  the  changes  are  caused  exclusively  by  the 
movement  of  the  coordinates  of  the  control  points  on  the 
boundaries  and  the  form  of  the  final  shape  is  pre-decided  in 
the  selection  of  the  initial  configuration.  If  a new  design 
is  desired,  it  would  be  necessary  to  change  the  boundaries, 
the  control  points,  and  even  the  mathematical  formulation  and 
repeating  the  optimization  the  procedure.  In  the  rule-based 
approach,  the  only  user  intervention  is  adding  the  user- 
defined  constraint  elements  to  the  data  file.  The  process 
will  proceed  automatically  to  generate  different  designs. 
Constraints  and  Loads 

In  the  finite  element  method,  the  information  about 
constraints  and  loads  is  related  to  the  nodes  of  structure. 
To  prepare  data  for  the  constraints,  the  user  of  the  finite 
element  program  must  find  out  nodal  identifications  and 
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constrained  degree  of  freedom  at  the  nodes  corresponding  to 
the  constraints.  For  the  surface  loads,  the  user  has  to  find 
the  identification  for  the  nodes  that  carry  load  and  the 
equivalent  nodal  load  must  be  calculated  at  the  corresponding 
node.  The  point  load  in  the  finite  element  analysis  should 
acts  on  the  node  only.  If  the  point  load  acts  inside  element 
instead  of  on  the  node,  the  user  has  to  calculate  equivalent 
nodal  loads  for  all  element  nodes  caused  by  the  point  load. 
In  the  analysis  of  a simple  structure,  only  a short  period  of 
time  in  needed  to  accomplish  this.  However  for  a complex 
structure,  the  node  identification  is  more  difficult  to  find 
out  and  it  may  take  considerable  time  to  prepare  the  data  for 
the  finite  element  model.  This  inconvenience  may  be 
acceptable  if  the  analysis  is  to  be  performed  one  time  only. 

In  the  automated  shape  design,  any  one  of  these 
identifications  may  cause  difficulties  since  the  data  must  be 
prepared  again  after  every  design  iteration.  This  is  the 
reason  that  some  parameterized  boundary  approaches  of  shape 
optimization  use  the  same  node  sequences  and  element 
identifications  in  every  design  iteration.  The  element  number 
and  node  sequence  are  constrained  without  any  changes  and  only 
the  coordinates  of  nodes  can  be  changed  to  move  the  control 
points  on  the  boundary.  The  drawback  of  this  is  that  the  ill- 
conditioned  finite  element  may  be  generated  during  the 
optimization  process,  because  the  movements  of  control  points 
depends  on  the  optimization  search  without  consideration  of 
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the  condition  of  the  finite  elements.  If  this  condition  is  to 
be  avoided,  the  remesh  algorithm  must  be  applied  to  the  system 
and  the  control  point  must  be  exactly  located  on  the  node.  The 
identification  of  control  points  also  should  be  updated  in 
every  design  iteration. 

In  this  study,  the  system  employs  the  mesh  generator  to 
remesh  the  design  domain  for  every  design  iteration.  The 
element  identifications  and  node  sequences  are  different  in 
every  design  iteration.  The  constraint  elements,  surface 
loading  elements  and  point  load  elements  are  included  in  the 
program  to  update  the  boundary  conditions  and  loads  for  the 
finite  element  model  during  the  optimization  process.  This  is 
previously  discussed  before  in  the  chapter  three.  The 
boundary  representation  model  is  used  in  the  program  to 
describe  the  shape  of  structure.  The  constraint  and  surface 
loading  elements,  like  the  boundary  of  the  structure,  use 
geometric  descriptions  combined  with  physical  data  to 
represent  the  physical  model  of  boundary  conditions  and  loads. 
The  user  does  not  need  to  find  out  the  node  identifications  or 
calculate  the  equivalent  nodal  forces  for  the  finite  element 
analysis.  The  user  has  to  describe  only  the  boundary 
conditions  and  loads  with  the  following  characteristics: 

(1)  Geometric  properties.  The  constraints  and  surface 
loads  may  be  represented  by  a line  or  arc.  For  the 
constraints  this  represents  where  the  boundary  is 
constrained  in  nodal  degree  of  freedom.  For  surface 
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loads,  this  represents  the  location  on  the  boundary  where 
the  surface  load  acts.  A point  load  is  only  represented 
by  a node  within  the  design  domain. 

(2)  Physical  data;  The  data  for  the  constraints  are  the 
flags  to  indicate  which  nodal  degree  of  freedom  is 
constrained.  The  forces  of  a point  load  and  pressure  of 
surface  loading  in  every  nodal  degree  of  freedom  are  the 
necessary  physical  data  for  the  finite  element  model. 
The  program  will  prepare  all  the  necessary  data  for  the  finite 
element  calculations  automatically.  The  data  for  the 
constraints  include  the  constrained  node  identification  and 
degrees  of  freedom.  The  data  for  the  loads  are  the  node 
identifications  and  the  corresponding  nodal  forces  caused  by 
the  surface  loads  or  point  loads.  In  this  way,  the  entire 
design  process  becomes  a fully  automatic  procedure  with 
minimum  user  intervention. 

Finite  Element  Calculations 

The  finite  element  method  is  used  in  this  study  to 
discretize  the  structure  and  to  calculate  the  stress 
distribution.  The  results  of  the  finite  element  analysis, 
such  as  displacements,  stresses,  velocities,  and 
accelerations,  can  then  be  used  to  guide  the  optimization 
procedures.  The  structure  with  static  or  inertial  loads  the 
displacements  is  first  analyzed  by  the  finite  element  method. 
The  nodal  stresses  are  calculated  by  averaging  the  nodal 
stresses  of  the  same  node  which  are  used  by  different 
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elements.  The  distortion  energy  theory  is  applied  to 
calculate  the  element  stress.  The  element  and  the  nodal 
stresses  are  then  used  for  the  elimination  algorithm. 

For  the  dynamic  loading,  the  Wilson  theta  method,  an 
implicit  direct  integration  method,  is  used  to  solve  the 
governing  equations.  The  damping  force  is  not  included  in  the 
system  equation.  The  time  step  plays  an  important  role  in 
solving  the  dynamic  problem.  The  total  computation  time 
needed  is  dependent  on  the  accuracy  of  the  numerical 
integration  technique  as  well  as  the  number  of  system 
solutions  it  requires  per  integration  step.  In  each  time  step 
we  have  to  evaluate  the  acceleration,  velocity,  and 
displacement  vectors  for  the  finite  element  model.  There  can 
be  considerable  computational  time  needed  when  dealing  with  a 
large  system. 

In  order  to  obtain  an  adequate  solution,  a good  numerical 
technique  must  be  selected.  An  unconditionally  stable 
integration  scheme  is  the  first  requirement.  The  effect  of 
the  unconditionally  stable  integration  scheme  is  derived  from 
the  fact  that  in  order  to  obtain  accuracy  in  integration,  the 
time  step  can  be  selected  without  the  requirement  for  a time 
step  formula  of  critical  conditions.  The  Wilson  theta  method 
is  the  unconditionally  stable  integration  scheme  that  is 
selected  to  solve  the  dynamic  system  equation  in  this  study. 
The  next  important  factor  is  to  choose  the  time  interval.  As 
we  discussed  in  chapter  three,  a stability  criteria  can  be 
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written  as 


At 


1 1 
c^  l-2b 


(5.1) 


where  c is  the  speed  of  sound  in  the  medium  and  1 is  the 
smallest  dimension  of  an  element.  For  Wilson  theta  method,  b 
is  equal  to  1/3.  Using  a good  numerical  scheme  and  a proper 
time  interval  in  solving  dynamic  system  can  insure  adequate 
results  and  save  computational  time. 

In  the  wave  propagation  problem,  the  free  boundary  must 
maintain  a zero  stress  condition.  The  incident  wave  is 
reflected  at  the  free  boundary  with  the  same  value  as  the 
incident  wave,  but  the  stress  is  reversed  because  the 
direction  of  travel  is  reversed.  The  total  deflection  at  a 
free  end  is  therefore  doubled  by  the  superposition  of  the 
incident  and  reflected  waves,  while  the  two  stress  components 
cancel  each  other.  On  the  other  hand,  at  a fixed  boundary  the 
zero-displacement  condition  must  be  satisfied  for  the  incident 
and  reflection  waves.  In  this  condition  the  incident 
displacement  wave  and  reflected  displacement  wave  have 
opposite  signs,  but  the  incident  stress  wave  and  reflected 
stress  wave  have  the  same  signs.  The  two  deflections  at  a 
fixed  boundary  cancel  each  other  satisfying  the  required  zero- 
displacement  condition,  while  the  reflected  wave  produces  a 
doubling  of  stress.  The  continuous  change  of  bovindary  and 
shape  makes  the  wave  propagation  more  complex  after  each 
iteration.  The  element  stresses  vary  as  time  changes  from  a 
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low  stress  at  one  time  state  to  a very  high  stress  at  another 
time  state.  Since  the  peak  stress  is  the  factor  which 
determines  material  failure,  it  is  used  as  the  criterion  for 
choosing  the  candidate  elements  for  elimination  in  the  shape 
optimization  process  of  structures  subjected  to  impulsive 
loads . 

The  proper  length  of  the  response  time  necessary  to 
collect  enough  peak  stress  information  for  the  algorithm  has 
to  be  decided.  If  a short  time  period  is  selected,  the 
elimination  procedure  may  choose  the  wrong  candidate  to 
eliminate,  because  if  a longer  time  is  allowed,  the  element 
eliminated  may  have  a high  peak  stress  due  to  the  stress 
build-up  phenomenon  in  wave  propagation.  On  the  other  hand  an 
unnecessary  long  time  can  be  computational  costly. 

We  select  the  cycle  time  of  response  as  twice  the  time 
needed  for  the  stress  wave  to  travel  the  longest  distance 
between  two  nodes  in  the  structure.  We  use  four  cycle  times 
as  the  response  time  within  which  the  peak  stresses  for  all 
elements  are  expected  to  occur.  The  natural  damping  in  most 
systems  can  be  expected  to  cause  the  response  to  decay  and 
consequently  the  peak  stresses  in  the  elements  will  not  be 
expected  to  occur  after  four  response  cycles.  The  duration  of 
the  impulsive  load  acting  on  the  structure  is  selected  as  half 
of  the  cycle  time. 


119 


Rule-Based  Elimination 

After  we  calculate  the  element  stresses  from  the  finite 
element  analysis,  the  next  step  is  to  select  the  element  to  be 
removed  from  the  structure.  Many  investigators  use  control 
points  on  the  master  elements  or  the  parameterized  boundaries 
to  change  the  shape  [8,9] . The  number  of  master  elements  or 
parameterized  boundaries  can't  be  changed.  The  shape  can  be 
only  changed  by  changing  these  pre-defined  master  elements  or 
boundaries.  In  the  case  of  complex  loads  and  boundary 
conditions , it  may  not  be  easy  to  pre-determine  the  master 
elements  or  controlled  boundaries.  This  approach  also  limits 
the  possible  alternative  solutions  that  may  lead  to  a better 
shape.  The  advantage  of  the  elimination  approach  is  that 
there  is  no  restriction  on  the  shape  and  the  algorithm  can  be 
applied  to  generate  general  shape.  It  is  not  pre-selected  in 
the  initial  step  and  the  optimum  is  generated  iteratively  in 
a natural  way. 

Any  element  without  any  constraint  or  load  on  it  can  be 
eliminated  from  the  structure  and  we  call  it  a free  element. 
The  calculated  stresses  for  the  free  elements  construct  the 
histogram  for  the  elimination.  Using  the  criterion  discussed 
in  chapter  four,  we  can  determine  the  element  identifications 
and  these  elements  are  grouped  into  elimination  groups.  The 
elimination  groups  are  the  subdomain  of  structure  that  will  be 
cut  from  the  structure  in  the  search  for  the  optimum  shape. 
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Change  and  Update  of  Shape 

The  elimination  groups  and  the  shape  of  structure  are 
united  by  the  Boolean  set  operation.  This  controls  the 
removal  of  material  from  the  structure  in  the  form  of  the  data 
structure  of  computer  model.  There  are  many  conditions  which 
may  occur  in  the  process  of  grouping  and  Boolean  set 
operation,  such  as  point  connections  between  elimination 
groups  or  between  elimination  group  and  the  structure.  A 
general  algorithm  is  developed  to  handle  the  different  cases 
in  a two-dimensional  domain.  The  finite  element  is  the  basic 
vehicle  in  the  optimization  process. 

After  removal  of  elements  from  the  structure,  some 
boundaries  of  structure  may  become  zigzag.  The  Bezier  spline 
is  applied  to  smooth  the  boundary.  The  order  of  the  base 
function  is  limited  because  if  we  use  the  maximum  available 
order,  it  becomes  a Bezier  curve.  If  the  defining  polygon 
vertices  can  create  a curve  with  large  curvature,  the  Bezier 
curve  will  be  further  away  from  the  defining  polygon  vertices 
than  the  Bezier  spline  with  the  proper  order  of  base  function. 
The  defining  polygon  vertices  used  in  this  study  are  the  nodes 
on  the  boundary  which  has  no  constraints  or  loads.  If  we  use 
a Bezier  curve  the  distance  between  the  control  polygon  and 
the  curve  is  generally  quite  large.  This  means  that  a larger 
area  between  the  control  polygon  and  the  curve  will  be  added 
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convergence  process  by  the  elimination  algorithm.  However  the 
Bezier  spline  is  flexible  to  a control  curved  shape  and  it  is 
chosen  to  be  used  in  the  program. 

Ending  the  Program 

A Criterion  should  be  given  to  end  the  optimization 
process.  In  theory  there  is  a point  where  further  elimination 
would  be  undesirable.  The  conditions  used  to  stop  the 
execution  of  the  algorithm  are  as  follows; 

( 1 ) The  average  stress  is  greater  than  seventy 
percentage  of  the  allowable  stress. 

(2)  No  further  elimination  is  possible  in  the  shape 
optimization  process. 

(3)  The  relative  changes  in  the  objective  function  is 
less  than  0.025  for  three  consecutive  design  iterations. 

(4)  The  maximum  number  of  iterations  is  violated. 

The  flow  chart  of  the  program  is  illustrated  in  figure  5-1. 
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Figure  5-1  Flow  chart  of  program  procedures 


CHAPTER  6 

APPLICATION  TO  STATIC  PROBLEMS 
Introduction 

In  order  to  illustrate  two-dimensional  applications  of 
the  shape  optimization  algorithm  by  elimination,  several  cases 
of  static,  inertial,  and  impact  loading  are  studied.  The 
initial  stress  distribution,  intermediate  shapes,  final 
result,  and  volume  changes  are  presented  for  the  treated 
cases.  Four  example  of  shape  optimization  under  static  loading 
are  discussed  in  this  chapter.  They  are  as  follows: 

(1)  Simple  structure  with  triangular  "load-support” 
configuration . 

(2)  Cantilever  structure. 

(3)  Statically  determinate  structure  with  multiple  loads. 

(4)  Statically  indeterminate  structure. 

All  cases  in  this  study  use  the  same  material.  The  material 
used  is  steel  with  the  following  properties: 

- Young's  modulus  30  000  000  psi 

- Poisson's  ratio  0.3 

- allowable  stress  40  000  psi 

The  allowable  stress  is  used  as  the  upper  constraint.  The 
structure  is  designed  such  that  no  part  of  it  exceeds  this 
stress.  The  surface  loading  in  the  treated  cases  is  set  equal 
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to  the  maximum  allowable  stress  in  most  cases.  All  cases  are 
started  with  a simple  initial  configuration,  a primitive 
shape,  which  satisfies  the  boundary  and  load  conditions  of  the 
design.  The  goal  of  the  optimization  program  is  to  reduce  the 
weight  of  the  structure  and  to  approach  as  closely  as  possible 
a uniform  stress  distribution. 

Case  1 - Simple  Triangular  Structure 

The  triangular  structure  is  the  simplest  structure  and  is 
commonly  used  as  the  basic  unit  to  construct  complex  trusses. 
The  initial  data  for  this 
case  is  shown  in  table  6-1. 

In  order  to  describe  the 
boundary  conditions  for  the 
finite  element  model,  the 
simply  support  is  modeled  by 
two  small  line  constraint 
elements,  one  with  constrained  degree  of  freedom  in  both  x and 
y directions  and  the  other  has  one  constrained  y degree  of 
freedom.  The  surface  load  is  model  by  a line  pressure  element 
with  pressure  equal  to  -40000  psi  in  y direction.  The  initial 
configuration  of  the  shape  used  for  starting  the  design 
procedure  is  rectangular  and  the  boundary  conditions  are  shown 
in  figure  6-1. 

The  stress  contours  for  the  initial  shape  are  plotted  in 
figure  6-2.  From  the  distribution  of  the  stress  contours,  we 
may  be  able  to  speculate  about  the  final  shape  of  the 


Table  6-1  Data  of  case  1 


length 

2 

width 

2 

meshes 

20  X 20 

load 

40,000  psi 
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optimized  structure.  The  stress  contours  determined  by  the 
finite  element  calculation  are  shown  in  figure  6-3  for  the 
third  iteration  and  the  corresponding  elements  which  are 
candidate  for  elimination  by  optimization  algorithm  are  shown 
in  figure  6-4.  During  the  process  of  optimization,  the 
eliminated  elements  are  always  located  in  the  area  of  the  low 
stress  contours  as  would  be  expected.  This  is  the  basic 
characteristic  of  elimination  concept.  Since  it  is  not  easy 
to  determine  the  elements  corresponding  to  the  stress 
contours,  the  use  of  the  histogram  in  this  study  embadies  the 
essence  of  the  stress  contours.  The  stress  contours  give  a 
display  of  the  stress  gradients  and  the  histogram  in  this 
study  is  constructed  by  dividing  the  elements  into  different 
groups  by  stress  levels  in  declining  orders  which  gives  a 
representation  of  the  stress  gradients  within  the  structure. 
Using  the  histogram,  we  can  obtain  more  statistical 
information  about  the  stress  distribution.  We  always 
eliminate  the  lowest  stressed  elements  and  insure  that  the 
boundary  constraints  are  not  violated.  Although  we  do  not 
eliminate  by  following  the  stress  contours,  the  outline  of  the 
eliminated  material  generally  follows  the  desired  low  stress 
contours.  The  eliminated  elements  which  are  identified  by  the 
shade  areas  in  figure  6-4,  are  removed  from  the  structure  by 
using  the  union  Boolean  set  operation.  After  the  union 
operation,  the  outline  of  the  structure  is  smoothed  by  the 
Bezier  spline  and  the  result  is  shown  in  figure  6-5.  This 
smooth  shape  is  used  for  the  next  optimization  iteration. 
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The  final  based  on  the  criterion  used  to  end  the  process 
shape  is  shown  in  figure  6-6  and  figure  6-7.  We  can  see  the 
stress  distribution  is  guite  uniform.  The  average  stress 
history  is  displayed  in  figure  6-8.  The  average  stress 
changes  from  10004  psi  to  23868  psi  and  the  maximum  stress  in 
the  free  elements  is  35710  psi  which  is  less  than  the 
allowable  stress  of  40000  psi. 

The  history  of  the  volume  of  the  structure  for  the 
different  iterations  is  shown  in  figure  6-9.  The  volume  of 
the  structure  changes  from  4.0  to  1.55  which  corresponds  to  a 
61.25  percent  of  reduction  in  volume.  This  simple  model  of 
structure  and  loads  illustrate  the  basic  concept  and  the 
feasibility  of  structural  optimization  by  elimination. 

In  order  to  insure  that  the  element  stress  does  not 
exceed  the  allowable  stress,  the  criterion  stress  of  ending 
the  search  is  seventy  percent  of  allowable  stress.  If  we  use 
a finer  mesh  in  the  finite  element  model,  the  criterion  stress 
can  be  increased  and  we  can  consequently  eliminate  more 
material  without  violating  the  allowable  stress. 

Case  2 - Cantilever  Structure 
The  initial  shape  and  mesh  in  this  case  are  the  same  as 
in  the  previous  case.  The  boundary  conditions  are  changed  to 
satisfy  the  constraint  condition  of  the  specific  cantilever 
structure.  The  degree  of  freedom  in  the  x and  y directions 
are  constrained  accordingly.  The  load  is  represented  by  a 
normal  pressure  acting  on  the  upper  right  corner  of  structure. 
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The  initial  configuration  and  the  boundary  conditions  are 
shown  in  figure  6-10.  The  initial  stress  distribution  is 
shown  in  figure  6-11.  During  the  optimization  process,  the 
elements  eliminated  from  the  structure  are  located  in  the  low 
stress  region  and  the  change  of  the  shape  generally  follows 
the  low  stress  contours. 

It  should  be  noted  that  when  we  eliminate  elements  with 
low  stress  level,  the  elements  carrying  high  stress  always 
stay  the  same  at  the  initial  steps  for  the  structure  which  is 
overdesigned  while  some  regions  of  structure  carry  low 
stresses.  The  elimination  speed  may  not  be  fast  enough  to 
remove  the  low  stress  regions  and  stress  distribution  may  not 
changed  significantly  in  the  high  stress  regions.  We  can  see 
this  in  the  change  from  first  iteration  shown  in  figure  6-11, 
to  iteration  three  and  four  shown  in  figure  6-12  and  figure  6- 
13  respectively. 

When  these  region  are  gone,  the  elimination  begins  to 
change  the  shape  and  the  stress  distribution.  Some  low  stress 
regions  may  begin  to  appear  in  the  region  that  is  previously 
loaded  with  higher  stress,  as  shown  in  figure  6-14.  The 
stress  distribution  in  some  of  the  high  stress  areas  may  also 
increase  as  shown  in  figure  6-15.  The  elimination  procedure 
approaches  the  final  result  gradually. 

The  final  shape  in  this  case  is  shown  in  figure  6-16  and 
figure  6-17.  We  can  see  that  the  final  result  is  influenced 
by  the  high  stress  distribution  regions  in  the  previous 
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iterations  such  as  that  shown  in  figure  6-15.  We  may  see  that 
the  final  structure  rises  two  triangles  to  carry  the  load 
which  would  be  case  when  using  a plane  truss.  If  we  want  to 
use  this  result  to  construct  a plane  truss  with  pin  joints  we 
may  use  two  triangular  truss  structures  based  on  the  final 
configuration  of  the  optimization  process. 

The  average  stress  changes  from  9693  psi  to  23545  psi. 
The  volume  of  the  structure  is  reduced  by  54.32  percent.  The 
volume  of  the  structure  may  be  reduced  more  by  using  a finer 
mesh  for  the  finite  element  calculations.  When  the 
optimization  iterations  approach  the  final  result,  the  stress 
becomes  higher  and  more  uniformly  distributed  over  the 
structure.  At  this  point  the  elimination  procedure  must  be 
slowed  down  because  there  are  only  few  elements  available  for 
elimination  which  will  not  affect  the  smoothed  boundary.  A 
finer  mesh  can  help  in  allowing  the  process  to  continue 
further. 

Case  3 - Multiple  Load  structure 

For  testing  the 
algorithm  in  a more  complex 
structure,  three  loads  are 
applied  on  the  bottom  of 
structure.  The  design  data 
is  given  in  table  6-2  and  the 
initial  configuration  and  stress  distribution  are  shown  in 
figure  6-18  and  figure  6-19  respectively. 


Table  6-2  Data  for  case  3 


length 

4 

width 

2 

meshes 

32  X 16 

loads 

40  000  psi 
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Plane  stress  analysis  using  the  finite  element  method  is 
utilized  in  this  case  as  in  all  the  considered  cases.  The 
stress  distribution  in  the  initial  configuration  is  shown  in 
figure  6-19.  The  final  optimum  shape  and  stress  distribution, 
as  shown  in  figure  6-20  and  figure  6-21,  can  not  be  used  to 
represent  a truss  structure  since  the  structure  can  carry 
bending  load  while  the  truss  element  can  only  carry 
longitudinal  forces.  The  result  of  this  case  can  not  be 
directly  used  as  a truss  structure.  It  may  be  used  however  as 
a guide  for  developing  the  configuration  of  the  truss  as  shown 
in  figure  6-22. 

The  average  stress  changes  from  9743  psi  to  21640  psi  and 
the  maximum  element  stress  is  35191  psi.  The  volume  of  the 
structure  is  reduced  by  56.41  percent.  The  volume  can  be 
further  reduced  if  a finer  mesh  is  used  to  slow  down  and 
refine  the  elimination  process. 

Case  4 - Statically  Indeterminate  Structure 

The  algorithm  of  shape  optimization  by  elimination  and 
refinement  is  a more  general  approach  than  the  mathematical 
programming  methods  since  the  later  have  to  undergo  a change 
in  formulation  each  time  the  problem  changes.  This  is 
illustrated  in  this  case  by  the  problem  of  optimizing  a 
statically  indeterminate  structure  with  boundary  conditions  as 
shown  in  figure  6-23.  The  stress  distribution  for  initial 
configuration  in  this  case  is  shown  in  figure  6-24. 
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The  finite  element  model  and  elimination  algorithm  be 
readily  used  to  optimize  the  structure.  There  is  no  other 
formulation  needed  when  we  change  from  one  problem  to  another. 
This  is  a definite  advantage  of  elimination  algorithm  and  it 
makes  it  ideally  suited  for  use  in  expert  system  for 
structural  optimization  in  the  future.  All  that  is  needed  to 
find  the  optimum  shape  of  structure  is  to  describe  the  loads, 
support  conditions,  geometric  constraints  and  the  finite 
element  model  for  a primitive  geometry  can  be  automatically 
constructed  and  the  algorithm  can  then  be  used  to  generate  the 
optimum  shape.  The  final  shape  and  stress  distribution  for  the 
considered  case  are  shown  in  figure  6-25  and  figure  6-26 
respectively. 

The  finite  element  modeling  used  in  the  shape 
optimization  algorithm  does  not  require  any  formulation  when 
dealing  with  a complex  structure.  It  is  therefore  possible  to 
solve  any  general  shape  optimization  problem  by  the 
combination  of  a rule-based  elimination  algorithm  and  a finite 
element  method  with  automatic  mesh  generation  as  illustrated 
in  the  previous  cases.  If  we  utilize  mathematical  programming 
methods  in  shape  optimization,  the  problem  formulation  must  be 
changed  to  fit  the  case.  The  generality  of  this  method  makes 
it  a powerful  tool  for  solving  complex  problem  without  the 
need  for  specialized  or  sophisticated  formulation. 
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Plan^^tr— Sttnjctur*  Shapa  — Iteration  Nunbr  t 1. 

Figure  6-1  Initial  shape  of  simple  supported 
structure  under  static  load 


averaae  stress  :l6004 .334* 


Figure  6-2  Initial  stress  distribution  of  simple 
supported  structure  under  static  load 
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Figure  6-3  Stress  distribution  of  iteration  three 
of  case  one 


case  one 
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Plana_Str«ss_CAs«l  Bourtdaru  Loop  — Iteration  Nu»«ber  : 4 

Figure  6-5  New  boundary  for  iteration  four  of  case 
one 
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Figure  6-7  Stress  distribution  of  final  shape  of 


case  one 
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static  load 


Figure  6-11  Stress  distribution  of  initial  shape 
of  case  two 
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Figure  6-12  Stress  distribution  of  iteration  3 of 
case  two 


Figure  6-13  Stress  distribution  of  iteration  four 
of  case  two 
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Figure  6-14  Stress  distribution  of  iteration  five 
of  case  two 


Figure  6-15  Stress  distribution  of  iteration  six 


of  case  two 
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Figure  6-17  Stress  distribution  of  final  shape  of 
case  two 
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Pl»n^_^tr»m«_Cam«3  Structure  Shap^  — ltT»tion  Nurib^r  ! 1 

Figure  6-18  Initial  shape  of  simple  supported 

structure  under  multiple  static  loads 


Figure  6-19  Stress  distribution  of  initial  shape 
of  case  three 
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Figure  6-21  Stress  distribution  of  final  shape  of 
case  three 
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construct  truss  structure 
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Plana_Strass_Pasa4  Structure  Shape  ~ Iteration  Nunber  : 1 

Figure  6-23  Initial  shape  of  statically 
indeterminate  case 


■auerafie  stress  : 10810.897 


Figure  6-24  Stress  distribution  of  initial  shape 
of  case  four 
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CHAPTER  7 

APPLICATION  TO  INERTIAL  LOAD  PROBLEM 
Introduction 

There  is  a host  of  practical  situations  in  which  the 
distribution  of  stress  manifests  symmetry  about  an  axis.  An 
example  of  this  class  of  problems  is  the  case  of  rotating  disk 
where  the  magnitude  and  distribution  of  the  stress  is  an 
important  factor  in  the  design  of  high  speed  machinery.  A 
rotating  disk  subjected  to  inertial  force  is  considered  to 
illustrate  the  variety  of  problems  where  structural 
optimization  can  be  accomplished  by  the  elimination  algorithm. 
When  optimizing  structures  subjected  to  static  loads  as  in  the 
previous  chapter,  we  can  treat  any  structural  optimization 
problem  by  using  a appropriate  initial  finite  element  model 
and  use  the  elimination  algorithm  without  further  mathematical 
formulation.  It  is  a very  convenient  way  to  solve  different 
problem.  This  is  also  the  case  when  we  change  the 

optimization  problem  from  the  plane  stress  case  to  the 
axisymmetric  case  and  the  only  thing  we  have  to  do  is  to 
construct  the  finite  element  model  for  the  axisymmetric  case. 
The  optimal  design  of  a rotating  disk  is  studied  as  an  example 
of  axisymmetric  problems.  The  energy  of  distortion  theory  is 
still  used  for  the  design  of  the  axisymmetric  cases. 
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Case  5 - Comparison  with  Hyperbolic  Profile 
The  maximum  stress  in  a rotating  disk  with  uniform 
thickness  is  known  to  occur  at  the  innermost  radius.  This 
explains  the  general  shape  of  many  disks:  thick  near  the  hub 
and  taping  down  in  thickness  toward  the  outermost  rim.  This 
not  only  has  an  effect  of  reducing  weight,  but  also  results  in 
lower  rotational  inertia.  Assuming  that  the  profile  at  any 
radial  section  is  represented  by  the  general  hyperbolic 
function 

t=t^r'®  (7.1) 

where  t^  is  a constant  and  s is  a positive  number,  the  shape 
of  curve  clearly  depends  on  the  value  selected  for  s.  For 
example,  for  s=l,  the  profile  is  that  of  an  equilateral 
hyperbolic.  The  constant  t^  represents  the  thickness  at  a 
radius  equal  to  unity.  If  the  thickness  at  r=a  is  t^  and 
thickness  at  r=b  is  t^^,  the  hyperbolic  curve  is  a curve  can  be 
evaluated  from 


to 


tia ^ /b\s 

t^b-"  la) 


(7.2) 


and  solving  for  s.  In  this  case,  we  use  thickness  tj=3  at  r=l 
and  thickness  t^j=0 . 5 at  r=6 , then 


h = -J- 

to  0.5 


(7.3) 


and  s can  be  solved  as  s=l. 
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In  order  to  compare  the  result  from  the  optimization 
algorithm  with  the  hyperbolic  curve  we  use  the  example  where 
the  thickness  t-=3  at  r=l  and  the  thickness  t^=0.5  at  r=6.  In 
this  case  s is  equal  to  one.  The  optimization  problem  is 
started  in  this  case  with  a rectangular  shape  at  the  initial 
configuration.  The  thickness  is  equal  to  three  where  the 
uniform  initial  shape  is  shown  in  figure  7-1  and  its  stress 
distribution  is  shown  in  figure  7-2.  We  can  see  that  the 
stress  increases  as  the  radius  decreases. 

The  elimination  algorithm  for  structural  optimization  is 
employed  to  determine  the  optimum  shape.  The  final  shape  is 
shown  in  figure  7-3.  The  stress  distribution,  as  shown  in 
figure  7-4,  is  much  more  uniform  than  in  the  previous  case  . 

In  order  to  compare  the  hyperbolic  profile  with  the 
optimization  result,  the  profile  of  optimization  result  is 
superimposed  in  the  same  graph  on  the  hyperbolic  profile  as 
shown  in  figure  7-5.  It  is  interesting  to  note  that  the  two 
profiles  are  very  close  to  each  other. 

The  above  case  is  repeated  with  a constrained  hub  with 
inner  radius  of  one  inch  and  outer  radius  of  two  inches.  The 
initial  configuration  and  final  results  for  this  case  are 
given  in  figure  7-6  and  figure  7-7.  The  hyperbolic  profile  in 
this  case  is  superimposed  on  the  optimal  profile  as  shown  in 
figure  7-8. 
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Case  6 - Optimal  Design  with  Different  Configurations 

The  rotating  disk  with  given  hub  and  rim  configurations 
is  considered  in  this  case.  The  size  of  the  hub  and  rim  are 
constraints  that  can  not  be  changed.  These  are  rectangular 
regions  and  this  kind  of  constraint  can  be  defined  by  a 
constraint  element  of  a rectangular  box  in  the  program.  The 
purpose  of  this  constraint  is  to  keep  the  material  in  that 
region  without  being  candidate  for  elimination.  There  is  no 
constrained  degrees  of  freedom  anywhere  inside  the  box.  The 
flags  used  to  indicate  the  constrained  degree  of  freedom  are 
set  to  zeros.  The  meaning  of  this  to  the  program  is 
equivalent  to  the  rule  that  the  element  inside  box  will  be 
kept  and  any  degree  of  freedom  inside  the  box  is  free  to  move. 

The  rotation  speed  in  all  the  cases  considered  in  this 
chapter  is  10000  rpm.  Material  used  is  the  same  as  material 
used  in  the  static  cases.  We  start  the  optimization  process 
with  a flat  rotating  disk.  In  order  to  keep  material  along 
the  center  line  at  the  junction  with  the  rim  from  being 
eliminated  a line  constraint  element  is  defined  and  the  reason 
for  this  constraint  element  will  be  discussed  later.  The 
initial  shape  and  its  stress  distribution  are  shown  in  figure 
7—9  and  figure  7—10  for  the  case  when  the  hub  and  rim  are 
given  as  rings  with  thickness  equal  to  3 inches. 

The  final  results  for  the  optimal  shape  are  shown  in 
figure  7-11  and  figure  7-12.  At  the  initial  step  the  stress 
distribution  is  shown  to  increase  as  the  radius  decreases. 
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The  stress  distribution  in  the  connecting  section  between  hub 
and  rim  in  the  optimum  shape  is  quite  uniform  as  shown.  The 
stress  values  the  rim  are  low  as  a result  of  the  imposed 
constraint  keep  all  the  material  there  no  matter  how  low  the 
stress  is  on  any  of  the  elements. 

Seireg  [26]  has  investigated  the  optimal  configurations 
with  different  rim  and  hub  constraints  by  mathematical 
programming  optimization.  The  result  for  the  same  conditions 
as  given  in  [26]  is  shown  in  figure  7-13  and  figure  7-14 
respectively . 

In  order  to  compare  the  two  results,  the  optimum  profiles 
determined  by  both  methods  are  superimposed  in  the  same  plot 
as  shown  in  figure  7-15.  We  can  see  that  the  two  results  are 
very  close.  The  profile  created  by  elimination  algorithm  is 
not  smooth.  We  can  use  a finer  mesh  to  achieve  a smooth 
boundary  or  we  can  smooth  the  boundary  by  curve  fitting  in 
according  to  the  output  from  computer. 

In  this  case  only  one  constraint  element  is  used  in  the 
center  of  the  connecting  section  between  the  web  and  the  rim. 
The  reason  for  this  is  to  compare  the  elimination  algorithm 
with  other  techniques  assuming  a single  web.  Another  reason 
is  to  demonstrate  the  flexibility  and  alternative  of  the 
elimination  concept  in  accommodating  the  designers 
specifications.  As  we  mention  in  chapter  four,  the 
parameterized  boundary  approach  can  only  change  the  position 
of  the  control  on  the  boundary  to  change  the  shape  and  it  does 
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not  create  a new  shape  for  the  boundary.  The  elimination 
concept  is  not  handicapped  by  this  constraint. 

The  same  initial  configuration  and  physical  model  as  in 
the  previous  case  is  now  considered  without  the  constraint 
element  in  the  center  of  the  connecting  section  between  the 
web  and  the  rim.  We  just  let  the  elimination  algorithm 
proceed  with  the  optimization  without  any  imposed  constraint 
on  the  structure.  The  optimum  configuration  in  this  case  is 
quite  different  from  the  previous  one  as  shown  in  figure  7-15 
and  figure  7-17.  Comparing  the  result  in  figure  7-11  and  7- 
17,  the  difference  in  the  total  volume  and  stress  is 
relatively  small  and  the  main  difference  is  the  actual 
configuration.  The  volume  in  the  hub  and  rim  can't  be 
changed.  The  volume  reduction  is  calculated  for  the 
connecting  section  between  the  hub  and  rim  only.  For  the 
first  optimal  design  shown  in  figure  7-11  the  volume  is 
reduced  by  60.94  percent  and  the  average  stress  is  13211  psi. 
The  maximum  stress  is  17305  psi  and  the  minimum  stress  is 
10926.  For  the  second  optimal  design  shown  in  figure  7-17, 
the  volume  is  reduced  by  74.77  percent.  The  average  stress  is 
17789  psi.  The  maximum  stress  is  20972  psi  and  the  minimum 
stress  is  14210  psi.  The  manufacturing  cost  for  the  first 
optimal  design  should  be  lower  than  that  for  the  second 
optimal  design.  This  example,  however,  illustrates  the 
potential  of  the  program  to  generate  options  for  the  user. 


151 


Another  case  is  considered  here  with  different  hub  and 
rim  configurations.  The  initial  configuration  and  the 
corresponding  stress  distribution  is  shown  in  figure  7-18  and 
figure  7-19.  The  final  results  of  the  optimization  are  shown 
in  figure  7-20  and  figure  7-21.  There  is  no  additional 
formulation  needed  for  this  case  since  all  the  necessary  rules 
and  analytical  procedures  are  the  same  for  all  cases. 

It  is  interesting  to  note  the  ring  which  is  generated 
immediately  after  the  hub.  This  is  similar  to  the  finding 
reported  by  Seireg  and  Surana  in  reference  [26].  If  this  ring 
is  trimmed  as  shown  in  figure  7-22  and  figure  7-23,  the 
maximum  stress  at  the  inner  radius  will  increase  from  17121 
psi  to  17151  psi  as  result. 


152 


Rotat  ir»gOsik_Hupar 

Structura 

Shapa 

I tarat ion 

Nunbar 

L 

Figure  7-1  Initial  shape  of  hyperbolic  profile 
case  without  hub 


Figure  7-2  Stress  distribution  of  initial  shape 
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Figure  7-3  Final  shape  of  hyperbolic  profile  case 


Figure  7-4  Stress  distribution  of  final  shape 
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Figure  7-6  Stress  distribution  of  initial  shape 


Figure  7-7  Stress  distribution  of  final  shape 
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Figure  7-9  Initial  shape  of  rotating  disc  with 
central  constraint 


Figure  7-10  Stress  distribution  of  initial  shape 
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Figure  7-11  Final  shape  of  rotating  disc  with 
central  constraint 


Figure  7-12  Stress  distribution  of  final  shape 
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Figure  7-13  Shape  created  by  mathematical 

programming  method  in  reference  [26] 


Figure  7-14  Stress  distribution  of  compared  shape 
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Figure  7-15  Comparison  between  the  results  by 
mathematical  programming  and 
elimination  algorithm 
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Figure  7-16  Final  shape  of  rotating  disc  without 
central  constraint 


Figure  7-17  Stress  distribution  of  final  shape 
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different  sizes  of  hub  and  rim 


Figure  7-19  Stress  distribution  of  initial  shape 
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different  sizes  of  hub  and  rim 


Figure  7-21  Stress  distribution  of  final  shape 
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Figure  7-22  Shape  of  rotating  disc  with  trimmed 
ring 
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Figure  7-23  Stress  distribution  of  rotating  disc 
with  trimmed  ring 


CHAPTER  8 

APPLICATION  TO  IMPACT  PROBLEM 
Introduction 

The  characteristics  of  the  impulsive  loading  is  that  the 
superposition  of  stress  waves  can  produce  much  larger  stress 
than  the  input  stress.  Finding  the  optimal  shape  under 
impulsive  loading  presents  a considerable  challenge.  The 
literature  devoted  to  the  general  problem  of  optimal  design  of 
structures  subjected  to  impact  loading  is  very  scare  if  not 
nonexistent.  In  this  chapter  we  will  use  the  same  elimination 
algorithm  to  optimize  two-dimensional  structures  subjected  to 
impact.  Unlike  the  static  or  inertial  loads  conditions,  the 
stress  wave  propagate  all  over  the  structure  for  a 
considerable  period  of  time.  In  this  period  of  time,  the 
stress  at  each  node  fluctuates  in  according  to  condition  of 
the  stress  wave  passing  through  node.  The  stress  history,  not 
the  stress  at  a specific  time,  is  the  controlling  design 
factor  in  the  optimization  process. 

Assuming  a continuous  two-dimensional  structure  subjected 
an  external  force  F(t)  applied  at  point  P,  this  external  force 
causes  a disturbance  inside  structure.  The  purpose  of  an 
analysis  is  to  compute  the  deformation  and  the  distribution  of 
stress  as  a function  of  the  spatial  coordinates  and  time.  If 
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the  largest  velocity  of  propagation  of  the  disturbance  is  c 
and  if  the  external  force  is  applied  at  t=0,  the  disturbed 
region  at  time  t=tl  is  generally  bounded  by  a circle  centered 
at  the  point  P with  radius  ctl.  Thus  the  entire  structure  is 
disturbed  at  time  t=r/c  where  r is  the  largest  distance  within 
the  structure  measured  from  point  P.  If  the  significant 
changes  in  F(t)  happen  over  a time  period  t^,  the  dynamic 
effect  are  important  if  t^  and  r/c  are  the  same  order  of 
magnitude.  For  the  excitation  sources  that  are  applied  for  a 
short  period,  the  effects  of  wave  motion  are  important  if  the 
time  interval  of  applications  is  of  the  same  order  of 
magnitude  as  the  characteristic  time  of  transmission  of  a 
disturbance  across  the  body.  In  the  cases  considered  in  this 
chapter,  an  external  force  F(t)  is  applied  to  the  structure  at 
point  P for  a period  t^=l/c  where  1 is  the  longest  distance 
between  point  P and  structure  boundary.  In  the  program  we 
define  a cycle  time  as  t^=2t^  that  is  twice  of  the  time  needed 
for  a stress  wave  to  travel  the  longest  distance  between  point 
P and  the  structure  boundary. 

An  important  and  useful  extension  of  Hooke's  law  is  the 
principle  of  superposition  which  enables  the  stresses  in 
interacting  waves  to  be  combined.  All  the  considered  cases 
are  the  elastic  structures  with  stresses  and  strains  linearly 
related.  The  principle  of  superposition  [59-62]  can  therefore 
be  applied  to  the  impact  problems  under  investigation. 
Superposition  frequently  takes  place  among  transient  stress 
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waves  and  leads  to  the  development  of  highly  localized 
concentrations  of  stress.  Interference  between  two  or  more 
elastic  waves  can  arise  within  a structure  in  several  ways. 
A common  situation  leading  to  interference  is  that  a wave 
strikes  two  or  more  free  surfaces  and  generates  new 
disturbance  at  each  surface  which  subsequently  interference. 

A specific  disturbance  combined  with  particular  boundary 
conditions  of  a structure  may  create  tension  or  compression 
waves.  When  two  stress  waves  meet  at  the  same  point  of 
structure,  one  in  tension  and  another  one  in  compression,  they 
will  annul  each  other  in  the  region  of  interaction.  This  will 
reduce  the  stress  in  that  region.  But  when  a pair  of 
compressional  waves  or  a pair  of  tensile  waves  meet  at  same 
point,  the  stress  will  add  up  in  the  region  of  interaction. 
This  transient  localized  stress  phenomenon  is  the  key  factor 
in  the  optimization  of  a structure  subjected  to  impact  loads. 
If  the  stress  of  a point  at  a specific  time  is  the  sum  of  more 
than  two  stress  waves  with  same  phase,  the  structure  will  fail 
when  the  stress  is  larger  than  the  allowable  value  according 
to  the  failure  criterion.  According  the  primary  task  in 
optimizing  a structure  subjected  to  impact  loads  is  to  reduce 
the  peak  transient  stresses  during  the  considered  period.  In 
the  developed  program,  we  evaluate  the  maximum  nodal  stresses 
for  the  entire  period  of  computation  and  this  is  the 
information  that  is  used  to  optimize  the  structure. 
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Optimizing  a Structure  Subjected  to  an  Axial  Impact 

The  first  case  is  a structure  subjected  to  a square  axial 
impulsive  load  with  a magnitude  of  20000  psi  acting  on  the 
structure  for  half  of  a cycle  time.  The  initial  configuration 
and  boundary  conditions  are  shown  in  figure  8-1.  The  stress 
distribution  of  the  initial  configuration  is  shown  in  figure 
8-2.  We  should  note  that  the  stress  distribution  displayed 
here  is  the  maximum  nodal  stresses  of  each  node  for  the  entire 
computational  period.  The  final  shape  and  the  corresponding 
stress  distribution  are  shown  in  figure  8-3  and  8-4 
respectively.  We  can  see  that  the  stress  distribution  in  the 
final  result  is  quite  constant  over  whole  structure. 

If  the  program  is  continued  with  the  objective  of 
reducing  the  volume  and  consequently  increasing  the  average 
stress,  the  solution  in  this  case  is  shown  in  figure  8-5  and 
8-6.  There  is  also  a small  notch  on  each  side  boundary  of  the 
final  shape.  This  is  created  by  removing  an  elements  on  those 
boundaries  and  the  low  order  Bezier  spline  is  used  for  smooth 
boundary.  Using  a higher  order  curve  to  get  more  smooth 
boundary  generating  a mesh  and  performing  a finite  element 
analysis  for  the  new  shape,  we  obtain  the  result  shown  in 
figure  8-7  and  figure  8-8.  The  results  compared  in  table  8-1. 

The  results  are  close,  but  we  can  see  that  the  original 
result  with  the  notch  is  better  than  the  result  from  the 
smoothed  shape  by  Bezier  curve.  The  original  results  also 
show  a smaller  difference  between  the  maximum  stress  and  the 
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Table  8-1  Comparison  among  different  shapes 


initial 

shape 

original 

result 

smoother 

shape 

maximum  element 
stress 

19226  psi 

26690  psi 

26644  psi 

minimum  element 
stress 

4579  psi 

17393  psi 

16187  psi 

mean  stress 

10937  psi 

19778  psi 

19634  psi 

stress  range 

14647  psi 

9297  psi 

10457  psi 

volume 

9.6000 
cubic  in . 

4.5612 
cubic  in. 

4.6068 
cubic  in. 

minimum  stress  and  the  volume  is  also  smaller.  A small  notch 
at  that  position  disturbs  the  wave  propagation  by  creating 
scatter  of  some  the  stress  waves.  The  results  indicate  the 
importance  of  the  notch  in  scattering  the  waves.  In  the 
design  of  structures  subjected  to  impact  properly  curved 
bovindary  may  reduce  the  peak  nodal  stress  by  avoiding  focused 
reinforcement  of  stress  wave. 

Optimizing  a Structure  Subjected  to  Lateral  Impact 

The  second  case  considered  for  dynamic  loading  is  the 
lateral  impact  on  a two-dimensional  structure.  The  initial 
configuration  and  the  stress  distribution  at  the  beginning  of 
the  optimization  process  are  shown  in  figure  8-9  and  8-10. 
The  final  shape  and  the  corresponding  stress  distribution  are 
shown  in  figure  8-11  and  8-12.  This  problem  represents  a more 
general  case  of  the  two-dimensional  analysis  and  the  final 
shape  is  also  more  complex.  In  analyzing  the  final  shape,  it 
is  noted  that  the  problem  is  very  sensitive  to  shape  changes 
in  the  final  steps  of  the  elimination.  A small  change  in  the 
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shape  may  result  in  local  focusing  of  stress  waves  and  may 
produce  significantly  higher  peak  stress.  In  the  case  of 
static  loads  or  body  forces  a small  change  in  shape  will  have 
only  a small  effect  on  the  stress  field  distribution,  which  is 
not  the  case  for  impact  loading. 

A summary  of  the  results  for  the  case  of  the  cantilever 
structure  subjected  to  impact  is  given  in  table  8-2. 


Table  8-2  Comparison  between  initial  shape  and  final  shape 


initial 

final 

maximum  element  stress 

30666  psi 

39351  psi 

minimum  element  stress 

8189  psi 

11445  psi 

mean  element  stress 

14374  psi 

20819  psi 

stress  range 

22547  psi 

27906  psi 

volume 

9.6  cubic  in. 

6.6  cubic  in. 
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Figure  8-2  Peak  stress  distribution  of  initial 
shape 
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runanirjg^l  Structure  Sh«p»  - Iteration  Nowbr  ; 20 

Figure  8-3  Final  shape  of  structure  subjected  to 
axial  impact 


Figure  8-4  Stress  distribution  of  final  shape 
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design  of  case  of  axial  impact 


Figure  8-6  Stress  distribution  of  shape 
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Figure  8-7  Smooth  shape  by  using  Bezier  curve 
instead  of  low  order  Bezier  spline 
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Figure  8-8  Stress  distribution  of  smooth  shape 
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Figure  8-10  Stress  distribution  of  initial  shape 
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Figure  8-12  Stress  distribution  of  final  shape 


CHAPTER  9 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  study  presented  in  this  thesis  investigates  the 
feasibility  of  using  a combination  of  a rule-based  elimination 
algorithm  and  the  finite  element  structural  modeling  in  a 
general  shape  optimization  procedure.  The  finite  element 
method  serves  as  a analytical  tool  of  analysis  for  this 
algorithm.  The  rule-based  of  elimination  and  mesh  generation 
procedures  proved  to  be  effective  in  all  the  cases  considered 
in  this  study.  The  cases  presented  in  chapters  six,  seven  and 
eight  are  intended  to  provide  different  classes  of  problems 
which  can  be  treated  by  this  approach. 

Based  on  the  previous  chapters , some  general  conclusions 
can  be  drawn  follows: 

(1)  The  algorithm  used  in  this  study  shows  excellent 
potential  for  providing  a general  tool  for  shape  optimization. 
Independent  of  the  type  of  load  and  boundary  conditions,  all 
that  is  needed  is  to  change  the  input  data  for  the  finite 
element  modeling  without  being  concerned  about  any  special 
formulation  of  the  problem.  This  is  particularly  helpful  in 
optimizing  the  shape  in  complex  problem  where  the  final 
configuration  cannot  be  predicted. 
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(2)  The  elimination  rules  are  very  important  in  the 
shape  optimization  process.  Different  criteria  of  eliminating 
material  from  the  structure  may  lead  to  different  results. 
Furthermore,  it  affects  the  total  iterations  needed  to 
complete  the  process.  The  criterion  for  elimination  is 
changed  based  on  to  the  stress  distribution  data  obtained  at 
each  step  during  the  optimization  process. 

(3)  We  always  expect  a zigzag  outline  for  the  structure 
after  the  elimination.  This  is  smoothed  by  the  application  of 
geometric  modeling  techniques.  After  changing  the  shape  in 
each  iteration,  the  new  outline  of  the  structure  is 
reconstructed  by  a combination  of  line,  arc,  and  Bezier  spline 
to  generate  a smooth  shape.  The  computer  modeling  technique 
for  the  two-dimensional  finite  element  mesh  generation  used  in 
this  thesis  is  quite  powerful  and  simple.  It  is  adapted  from 
geometric  modeling  of  two-dimensional  computer  graphics  and 
can  be  readily  extended  as  a general  technique  for  any  finite 
element  problem.  The  technique  helped  make  the  preparation  of 
data  for  finite  element  models  easy,  flexible,  and  fast. 
Furthermore,  the  physical  description  of  the  constraint  and 
load  elements  is  made  clear  practical  and  realistic. 

(4)  The  shape  generated  by  the  algorithm  in  the  static 
cases  is  directly  applicable  to  stamped  structures,  since  the 
elements  used  in  the  analysis  can  carry  shear  stress  and  this 
allows  the  structure  to  carry  bending  loads.  The  results  can 
form  the  basis  for  truss  design.  Since  truss  structures  can 


179 


only  carry  axial  load,  the  optimized  shape  can  provide  the 
outline  of  shape  for  the  layout  of  the  truss  structure. 

(5)  The  concepts  of  grouping  elements  for  elimination 
and  using  a Boolean  set  operation  make  the  elimination 
algorithm  free  of  any  restriction  on  changing  the  boundaries 
of  the  structure.  It  can  also  create  new  boundaries 
automatically.  This  can't  be  accomplished  by  the  mathematical 
programming  methods. 

(6)  In  the  rotating  disk  cases,  we  can  obtain  alternate 
solutions  depending  on  the  imposed  constraints  on  the 
elimination.  The  user  can  therefore  have  several  design 
options  to  select  from  based  on  considerations  of  cost  and 
case  of  manufacture. 

(7)  The  shape  of  structure  subjected  to  the  impact  loads 
is  quite  sensitive  to  the  boundaries  of  structure.  When  the 
design  proceeds  towards  the  optimal  shape,  the  elimination 
process  must  be  slowed  down.  This  also  applies  to  the  other 
cases,  but  it  is  more  important  for  the  dynamic  load  cases. 

(8)  The  split  line  criterion  of  mesh  generation  can  be 
improved  by  the  information  of  the  environment  about  a 
specific  point.  If  a general  split  criterion  is  defined,  the 
split  line  algorithm  may  automatically  generate  symmetry  in 
the  distribution  of  meshes  for  the  case  of  a symmetric  shape. 

In  order  to  construct  a sophisticated  system  which 
applies  the  concept  of  elimination  in  the  process  of 
structural  optimization,  we  must  integrate  geometric  modeling. 
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extensive  finite  element  functions,  and  a generalized 
criterion  of  elimination.  Some  recommendations  for  extending 
the  reported  study  are  given  in  the  following: 

(1)  The  modeling  stage  can  be  improved  by  incorporating 
the  geometric  modeling  technigues  used  in  three-dimensional 
CAD  systems  in  the  program  in  order  to  describe  the  shape  in 
a flexible  and  precise  way  for  all  problems.  Furthermore,  the 
advanced  techniques  of  finite  element  analysis,  such  as  error 
estimation,  adaptive  refinement,  should  be  included  in  the 
system.  This  integration  will  enable  the  algorithm  to  handle 
complex  structural  optimization  problem  more  easily  and 
precisely. 

(2)  Other  functions  of  the  finite  element  method  may  be 
included  to  solve  different  problems.  The  general  finite 
element  analysis  functions  from  one-dimensional  to  three- 
dimensional  problems  will  allow  the  system  the  to  treat  much 
more  complex  problems.  The  ability  of  the  finite  element 
method  to  solve  static  or  dynamic  problems  in  the  field  of 
elasticity,  plasticity,  and  thermal  stress  can  allow  us  to 
study  the  feasibility  of  the  elimination  concept  in  other 
field  of  applications. 

( 3 ) More  sophisticated  rules  and  criterion  for 
elimination  algorithm  can  be  developed.  The  rules  developed 
in  this  study  are  in  the  field  of  plane  stress  of  elastic 
structures,  rotating  disks,  and  wave  propagation.  The 
elimination  criterion  is  the  most  important  factor  in  the 
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optimization.  More  cases  should  be  studied  to  develop  a 
universal  general  criterion.  We  may  extend  the  concept  in  the 
field  of  eigenvalue  problems,  thermal  stress  analysis,  or 
plasticity  mechanics.  The  study  can  be  also  extended  to  deal 
with  problem  that  includes  more  than  two  field  of  loading  such 
as  combinations  of  static  forces  impulsive  loading  as  well  as 
steady  state  or  transient  thermal  stress. 

(4)  The  propagation  of  stress  waves  is  extremely 
difficult  to  predict  when  the  shape  of  the  structure  is 
changed.  Focusing  and  scattering  are  two  important  factors 
affecting  the  peak  stress  in  the  elements.  In  order  to  obtain 
uniform  peak  stress  distribution,  it  is  advisable  to  let  same 
stress  waves  pass  same  location  with  a small  time  delay.  Then 
the  superposition  of  total  stresses  becomes  high  enough 
without  violating  the  allowable  stress.  The  theory  about  how 
to  control  the  shape  to  reduce  the  peak  nodal  stress  should  be 
studied  and  for  use  in  the  elimination  procedure  to  improve 
the  optimal  design  of  structure  subjected  to  impact  loads. 

(5)  An  algorithm  for  adding  material  to  the  structure 
should  be  investigated  for  possible  improvement  of  the  process 
of  shape  synthesis  by  elimination  only.  We  may  use  the  result 
of  two  successive  iterations  as  the  basic  information  in  the 
addition  algorithm.  However,  it  is  possible  that  if  the 
design  criterion  is  improved,  there  is  no  need  to  add  material 
back  to  structure.  As  an  example  the  design  merit  is  worse  in 
the  current  iteration,  we  may  modify  the  shape  to  fall  between 
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the  shapes  of  these  two  iterations  and  restart  the  iteration 
from  the  new  modified  shape.  The  boundaries  in  region  with 
high  stress  value  may  be  shifted  by  a proper  distance  based  on 
a programmed  rule  to  add  material  back  to  structure. 

(6)  The  combination  of  the  elimination  group  concept  and 
the  Boolean  set  operation  provided  considerable  freedom  and 
flexibility  for  changing  the  shape.  However  there  are  a lot 
of  situations  to  deal  with  for  generalizing  the  process.  The 
algorithms  to  process  different  geometric  conditions  that  are 
created  by  the  elimination  optimization  can  be  further 
improved  and  extended  to  three-dimensional  problems. 


APPENDIX  A 

COMPUTER  RULES  FOR  SHAPE  OPTIMIZATION 
Find  Elimination  Candidate 

For  all  elements  in  the  design  domain,  the  procedures  listed 
below  are  followed  to  select  the  candidate  elements  for 
elimination: 

(1)  Evaluate  nodal  stresses  by  the  finite  element  method. 

(2)  Evaluate  element  stresses  using  distortation  energy 
theory . 

( 3 ) Identify  the  constrained  elements . 

( 4 ) Find  the  maximum  and  minimum  stresses  of  the 
unconstrained  elements . 

(5)  Define  the  group  number  at  this  step. 

(6)  Construct  the  histogram. 

(7)  Define  the  stress  distribution  by  mean  and  standard 
deviation. 

(8)  Find  the  maximum  increase  of  two  consecutive  groups 
with  stress  below  one  standard  deviation  of  mean. 

(9)  All  elements  below  that  are  candidates  for 
elimination. 

Grouping  Eliminated  Elements 

For  all  eliminated  elements,  the  number  of  adjacent  eliminated 
groups  connected  to  the  element  are  counted  and  the  following 
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rules  are  used  to  group  the  eliminated  elements  into 
eliminated  groups; 

(1)  If  there  is  no  adjacent  eliminated  group,  the 
element  is  used  to  form  a new  group. 

( 2 ) If  there  is  only  one  adj  acent  group , the  element  is 
joined  with  the  eliminated  group. 

( 3 ) If  there  are  more  than  one  adj  acent  groups , the 
element  is  jointed  with  one  of  the  connected  group,  then 
the  remaining  groups  are  joined  with  that  group. 

Determining  The  Boundary  of  The  Eliminated  Group 
Using  the  boundary  representation  model  for  the  eliminated 
groups,  the  boundary  of  the  eliminated  group  is  searched  and 
defined  as  follows: 

(1)  The  extreme  point  with  minimum  x coordinate  and 
minimum  y coordinate  is  selected  as  the  starting  point  to 
trace  the  boundary  of  the  group  and  it  is  defined  as  the 
first  bovindary  point. 

(2)  From  the  current  boundary  point,  the  number  of 
eliminated  elements  of  the  eliminated  group  sharing  the 
current  boundary  point  is  counted  and  the  next  boundary 
point  is  decided  according  to  the  following  rules: 

(a)  If  there  is  only  one  eliminated  element  using 
the  current  boundary  point,  the  next  boundary  point 
is  the  next  node  on  the  element  according  to  the 
nodal  seguence  of  the  eliminated  element. 

(b)  If  there  are  more  than  one  eliminated  element 
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sharing  the  current  boundary  point,  it  is 
necessary  to  find  the  single  edge  which  is  used  by 
one  eliminated  element  only  and  the  current 
boundary  point  is  one  of  the  end  points  of  the 
single  edge.  The  other  end  of  the  single  edge,  in 
counter-clockwise  sequence,  leads  to  the  next 
boundary  point. 

Preprocessing  For  Candidate  Groups  Before  Elimination 
In  order  to  simplify  the  conditions  occurring  in  Boolean 
operations,  the  eliminated  groups  are  processed  to  modify  some 
special  conditions  as  follows: 

( 1 ) The  point-connected  condition  is  checked  between  any 
two  eliminated  groups.  If  there  is  a point-connected 
between  two  eliminated  groups,  they  are  joined  into  one 
group . 

(2)  The  point-connected  condition  between  eliminated 
groups  and  shape  boundaries  is  also  checked.  If  there  is 
a point-connected  condition  between  any  eliminated  group 
and  shape  boiindary,  the  eliminated  group  is  modified. 

( 3 ) The  total  enclosure  is  also  checked  between  the 
eliminated  group  and  shape  boundary.  If  an  inside  hole 
is  totally  enclosed  by  the  eliminated  group,  the  outside 
boundary  of  eliminated  group  becomes  the  inside  hole  of 
the  structure. 

( 4 ) The  discontinuity  condition  is  checked  and  processed 
before  the  Boolean  operation.  If  there  is  a 


186 


discontinuity  condition,  the  following  rules  are  used: 

(a)  If  the  eliminated  group  connects  two 
constrained  region,  the  group  is  subdivided  into 
three  subgroups  and  the  central  group  is  kept  to 
maintain  connectivity  between  the  constrained 
regions . 

(b)  If  a disconnected  region  is  unconstrained  with 
a free  boundary,  it  is  removed  with  the  eliminated 
group . 

Update  Shape  Boundaries 

After  all  eliminated  groups  are  defined,  the  boundary  of  the 
shape  of  the  last  iteration  is  used  as  the  reference  shape  to 
update  the  shape.  The  following  rules  are  used  for  all 
eliminated  groups. 

(1)  The  intersection  between  the  eliminated  group  and 
the  shape  boundary  is  checked 

(2)  If  there  is  no  intersection  between  the  shape 
boundary  and  the  eliminated  group,  this  will  form  a new 
structural  boundary (a  hole). 

( 3 ) If  there  is  only  one  intersection  between  the  shape 
boundary  and  the  eliminated  group,  the  Boolean  operation 
is  applied  to  update  the  shape  boundary. 

( 4 ) If  there  are  more  than  one  intersection  between  two 
or  more  shape  boundaries  and  the  eliminated  group,  first 
the  Boolean  operation  is  applied  between  the  eliminated 
group  and  one  of  connected  shape  boundaries.  The  Boolean 
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operation  is  applied  again  between  this  new  updated  shape 
boundary  and  the  remaining  shape  boundaries. 

Generation  of  A New  Shape  After  Elimination 
The  boundary  of  the  shape  of  the  last  iteration  is  used  as  the 
reference  shape  to  generate  the  shape  after  applying  the 
elimination  algorithm.  The  following  rules  are  applied  to  the 
all  shape  boundaries. 

(1)  All  lines  and  arcs  on  the  boundaries  of  the  last 
iteration  are  located. 

(2)  All  points  on  the  newly  generated  boundary  after 
elimination  are  checked  to  determine  if  they  belong  to 
any  line  or  arc  of  the  reference  boundary  or  not.  If  the 
point  is  located  on  any  line  or  arc  of  the  reference 
boundary,  it  should  be  marked  with  the  identification  of 
the  corresponding  line  or  arc. 

(3)  The  points  with  line  or  arc  identification  will  be 
used  to  construct  a line  or  arc  on  the  new  shape 
boundary.  The  points  without  any  line  or  arc 
identification  are  the  points  which  are  previously 
located  on  the  Bezier  spline  of  the  reference  boundary  or 
are  the  newly  created  boundary  after  elimination.  All  of 
these  point  will  be  used  to  construct  a Bezier  spline  on 
the  new  shape  boundary. 

( 4 ) If  there  are  less  than  three  consecutive  points  with 
line  or  arc  identification  left  between  two  Bezier 
splines,  they  are  checked  for  constraints.  If  there  is 


188 


no  constraint  imposed  on  those  points,  they  can  be  set  to 
be  free  points  and  will  be  used  for  the  Bezier  spline 
instead  of  line  or  arc. 

(5)  The  common  point  between  two  consecutive  lines  or 
arcs  will  be  checked  for  special  conditions  and  the  line 
or  arc  identifications  at  the  common  point  are  rearranged 
such  that  the  program  can  trace  the  boundary  of  the  new 
shape . 

(6)  If  some  consecutive  points  have  the  same  line  or  arc 
identification,  they  will  be  used  to  construct  a line  or 
arc  for  the  new  shape.  The  consecutive  points  without 
any  line  or  arc  identification  will  be  used  to  construct 
the  Bezier  spline  for  the  new  shape. 


APPENDIX  B 

PROGRESSION  OF  SHAPE  MODIFICATION  FOR 
THE  CONSIDERED  EXAMPLES 

Progression  of  elimination  for  the  all  considered 
examples  is  illustrated  by  the  selected  sequence  of  results. 
The  stress  distributions  and  the  corresponding  histograms  for 
the  selected  iteration  number  are  given  in  figures  B-1  to  B-11 
respectively : 
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Iteration  9 


Iteration  14 


Figure  B-1  Simple  supported  structure  under  static  load 
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Iteration  3 


Iteration  4 


Iteration  10 


Iteration  18 


Figure  B-2  The  cantilever  beam  under  static  load 
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Figure  B-3  Simple  supported  structure  under  multiple  static 
load 
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Iteration  21 


Figure  B-4  Statically  indeterminate  structure  subjected  to 
static  load 
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Iteration  10 


Iteration  20  Iteration  30 


Iteration  46 


Iteration  58 


Figure  B-5  Rotating  disc  without  hub  constraint 
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Iteration  1 


Iteration  8 


Iteration  32 


Iteration  39 


Figure  B-6  Rotating  disc  with  hub  constraint 
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Iteration  1 


Iteration  15 


Iteration  25 


Figure  B-7  Rotating  disc  with  equal  thickness  rings  at  the 
hub  and  rim  with  elimination  constraint  on  center 
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Iteration  20 


Iteration  23 


Figure  B-8  Same  as  before  without  elimination  constraint  on 
center 
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Iteration  1 


Iteration  3 


Iteration  6 


Iteration  11 


Iteration  21 


Iteration  30 


Figure  B-9  Rotating  disc  with  constraint  on  the  hub  and  a rim 
with  large  width  than  the  hub 
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Iteration  4 


Iteration  10 


Iteration  1 


Iteration  7 


Iteration  15  Iteration  20 


Figure  B-10  Structure  subjected  to  axial  impact 
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Iteration  5 Iteration  9 


Iteration  11 
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Iteration  13 


Figure  B-11  Cantilever  structure  subjected  to  lateral  impact 
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